GLVis  v4.2
Accurate and flexible finite element visualization
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
vssolution3d.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-443271.
4 //
5 // This file is part of the GLVis visualization tool and library. For more
6 // information and source code availability see https://glvis.org.
7 //
8 // GLVis is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include <cstdlib>
13 #include <iostream>
14 #include <cmath>
15 #include <limits>
16 
17 #include "mfem.hpp"
18 #include "visual.hpp"
19 #include "palettes.hpp"
20 using namespace mfem;
21 using namespace std;
22 
23 
24 thread_local VisualizationSceneSolution3d *vssol3d;
25 extern thread_local GeometryRefiner GLVisGeometryRefiner;
26 
27 // Reference geometries with a cut in the middle, based on subdivision of
28 // GLVisGeometryRefiner in 3-4 quads. Updated when cut_lambda is updated, see
29 // keys Ctrl+F3/F4. We need these variables because the GLVisGeometryRefiner
30 // cashes its RefinedGeometry objects.
31 thread_local IntegrationRule cut_QuadPts;
32 thread_local Array<int> cut_QuadGeoms;
33 thread_local IntegrationRule cut_TriPts;
34 thread_local Array<int> cut_TriGeoms;
35 
36 // Definitions of some more keys
37 
38 std::string VisualizationSceneSolution3d::GetHelpString() const
39 {
40  std::stringstream os;
41  os << endl
42  << "+------------------------------------+" << endl
43  << "| Keys |" << endl
44  << "+------------------------------------+" << endl
45  << "| a - Displays/Hides the axes |" << endl
46  << "| A - Turns antialiasing on/off |" << endl
47  << "| c - Toggle colorbar and caption |" << endl
48  << "| C - Change the main plot caption |" << endl
49  << "| e - Displays/Hides the elements |" << endl
50  << "| E - Toggle the elements in the CP |" << endl
51  << "| f - Smooth/Flat/discont. shading |" << endl
52  << "| g - Toggle background |" << endl
53  << "| G - Export to glTF format |" << endl
54  << "| h - Displays help menu |" << endl
55  << "| i - Toggle cutting plane |" << endl
56  << "| I - Toggle cutting plane algorithm|" << endl
57  << "| j - Turn on/off perspective |" << endl
58  << "| k/K Adjust the transparency level |" << endl
59  << "| ,/< Adjust color transparency |" << endl
60  << "| l - Turns on/off the light |" << endl
61  << "| L - Toggle logarithmic scale |" << endl
62  << "| m - Displays/Hides the mesh |" << endl
63  << "| M - Toggle the mesh in the CP |" << endl
64  << "| o/O (De)refine elem, disc shading |" << endl
65  << "| p/P Cycle through color palettes |" << endl
66  << "| q - Quits |" << endl
67  << "| r - Reset the plot to 3D view |" << endl
68  << "| R - Reset the plot to 2D view |" << endl
69  << "| s - Turn on/off unit cube scaling |" << endl
70  << "| S - Take snapshot/Record a movie |" << endl
71  << "| t - Cycle materials and lights |" << endl
72  << "| u/U Move the level surface |" << endl
73  << "| v/V Add/Delete a level surface |" << endl
74  << "| w/W Move bdr elements up/down |" << endl
75  << "| x/X Rotate cutting plane (phi) |" << endl
76  << "| y/Y Rotate cutting plane (theta) |" << endl
77  << "| z/Z Translate cutting plane |" << endl
78  << "| \\ - Set light source position |" << endl
79  << "| Ctrl+o - Element ordering curve |" << endl
80  << "| Ctrl+p - Print to a PDF file |" << endl
81  << "+------------------------------------+" << endl
82  << "| Function keys |" << endl
83  << "+------------------------------------+" << endl
84  << "| F1 - X window info and keystrokes |" << endl
85  << "| F2 - Update colors, etc. |" << endl
86  << "| F3/F4 - Shrink/Zoom bdr elements |" << endl
87  << "| Ctrl+F3/F4 - Cut face bdr elements |" << endl
88  << "| F5 - Set level lines |" << endl
89  << "| F6 - Palette options |" << endl
90  << "| F7 - Manually set min/max value |" << endl
91  << "| F8 - List of subdomains to show |" << endl
92  << "| F9/F10 - Walk through subdomains |" << endl
93  << "| F11/F12 - Shrink/Zoom subdomains |" << endl
94  << "+------------------------------------+" << endl
95  << "| Keypad |" << endl
96  << "+------------------------------------+" << endl
97  << "| 1-9 Small rotation, reset with 5 |" << endl
98  << "| *,/ Scale up/down |" << endl
99  << "| +/- Change z-scaling |" << endl
100  << "| . - Start/stop spinning |" << endl
101  << "| 0/Enter - Spinning speed and dir. |" << endl
102  << "+------------------------------------+" << endl
103  << "| Mouse |" << endl
104  << "+------------------------------------+" << endl
105  << "| left btn - Rotation |" << endl
106  << "| middle btn - Translation |" << endl
107  << "| right btn - Scaling |" << endl
108  << "| left + Alt - Tilt |" << endl
109  << "| left + Shift - Spinning |" << endl
110  << "| right + Shift - Change light pos. |" << endl
111  << "| left + Ctrl - Spherical rotation |" << endl
112  << "| middle+ Ctrl - Object translation |" << endl
113  << "| right + Ctrl - Object scaling |" << endl
114  << "| left + Ctrl + Shift - z-Spinning |" << endl
115  << "+------------------------------------+" << endl;
116  return os.str();
117 }
118 
119 static void KeyiPressed()
120 {
121  vssol3d -> ToggleCuttingPlane();
122  SendExposeEvent();
123 }
124 
125 static void KeyIPressed()
126 {
127  vssol3d -> ToggleCPAlgorithm();
128  SendExposeEvent();
129 }
130 
131 void VisualizationSceneSolution3d::PrepareOrderingCurve()
132 {
133  bool color = draworder < 3;
134  order_buf.clear();
135  order_noarrow_buf.clear();
136  PrepareOrderingCurve1(order_buf, true, color);
137  PrepareOrderingCurve1(order_noarrow_buf, false, color);
138  updated_bufs.emplace_back(&order_buf);
139  updated_bufs.emplace_back(&order_noarrow_buf);
140 }
141 
142 
143 void VisualizationSceneSolution3d::PrepareOrderingCurve1(gl3::GlDrawable& buf,
144  bool arrows,
145  bool color)
146 {
147  gl3::GlBuilder builder = buf.createBuilder();
148 
149  // make the lines of the ordering curve thicker
150  double ThicknessFactor = 2.0;
151  double MS_Thickness = 2.0;
152  double LineWidth;
153  if (GetMultisample() > 0)
154  {
155  LineWidth = GetLineWidthMS();
156  SetLineWidthMS(MS_Thickness);
157  }
158  else
159  {
160  LineWidth = GetLineWidth();
161  SetLineWidth(ThicknessFactor*LineWidth);
162  }
163 
164  DenseMatrix pointmat;
165  Array<int> vertices;
166 
167  DenseMatrix pointmat1;
168  Array<int> vertices1;
169 
170  int ne = mesh->GetNE();
171  for (int k = 0; k < ne-1; k++)
172  {
173  mesh->GetPointMatrix (k, pointmat);
174  mesh->GetElementVertices (k, vertices);
175  mesh->GetPointMatrix (k+1, pointmat1);
176  mesh->GetElementVertices (k+1, vertices1);
177  int nv = vertices.Size();
178  int nv1 = vertices1.Size();
179 
180  ShrinkPoints(pointmat, k, 0, 0);
181  ShrinkPoints(pointmat1, k+1, 0, 0);
182 
183  double xs = 0.0;
184  double ys = 0.0;
185  double zs = 0.0;
186  for (int j = 0; j < nv; j++)
187  {
188  xs += pointmat(0,j);
189  ys += pointmat(1,j);
190  zs += pointmat(2,j);
191  }
192  xs /= nv;
193  ys /= nv;
194  zs /= nv;
195 
196  double xs1 = 0.0;
197  double ys1 = 0.0;
198  double zs1 = 0.0;
199  for (int j = 0; j < nv1; j++)
200  {
201  xs1 += pointmat1(0,j);
202  ys1 += pointmat1(1,j);
203  zs1 += pointmat1(2,j);
204  }
205  xs1 /= nv1;
206  ys1 /= nv1;
207  zs1 /= nv1;
208 
209  double dx = xs1-xs;
210  double dy = ys1-ys;
211  double dz = zs1-zs;
212  double ds = sqrt(dx*dx+dy*dy+dz*dz);
213 
214  if (color)
215  {
216  double cval = minv+double(k)/ne*(maxv-minv);
217  MySetColor(builder, cval, minv, maxv);
218  }
219 
220  if (arrows)
221  {
222  Arrow3(builder,
223  xs,ys,zs,
224  dx,dy,dz,
225  ds,0.05);
226  }
227  else
228  {
229  Arrow3(builder,
230  xs,ys,zs,
231  dx,dy,dz,
232  ds,0.0);
233  }
234  }
235 
236  if (GetMultisample() > 0)
237  {
238  SetLineWidthMS(LineWidth);
239  }
240  else
241  {
242  SetLineWidth(LineWidth);
243  }
244 }
245 
246 void VisualizationSceneSolution3d::CPPrepare()
247 {
248  PrepareCuttingPlane();
249  PrepareCuttingPlaneLines();
250 }
251 
252 void VisualizationSceneSolution3d::CPMoved()
253 {
254  CPPrepare();
255  if (cplane == 2)
256  {
257  Prepare();
258  PrepareLines();
259  }
260 }
261 
262 static void KeyxPressed()
263 {
264  vssol3d -> CuttingPlane -> IncreasePhi();
265  vssol3d -> FindNodePos();
266  vssol3d -> CPMoved();
267  SendExposeEvent();
268 }
269 
270 static void KeyXPressed()
271 {
272  vssol3d -> CuttingPlane -> DecreasePhi();
273  vssol3d -> FindNodePos();
274  vssol3d -> CPMoved();
275  SendExposeEvent();
276 }
277 
278 static void KeyyPressed()
279 {
280  vssol3d -> CuttingPlane -> IncreaseTheta();
281  vssol3d -> FindNodePos();
282  vssol3d -> CPMoved();
283  SendExposeEvent();
284 }
285 
286 static void KeyYPressed()
287 {
288  vssol3d -> CuttingPlane -> DecreaseTheta();
289  vssol3d -> FindNodePos();
290  vssol3d -> CPMoved();
291  SendExposeEvent();
292 }
293 
294 static void KeyzPressed()
295 {
296  vssol3d -> CuttingPlane -> IncreaseDistance();
297  vssol3d -> FindNodePos();
298  vssol3d -> CPMoved();
299  SendExposeEvent();
300 }
301 
302 static void KeyZPressed()
303 {
304  vssol3d -> CuttingPlane -> DecreaseDistance();
305  vssol3d -> FindNodePos();
306  vssol3d -> CPMoved();
307  SendExposeEvent();
308 }
309 
310 static void KeymPressed()
311 {
312  vssol3d -> ToggleDrawMesh();
313  SendExposeEvent();
314 }
315 
316 static void KeyePressed()
317 {
318  vssol3d -> ToggleDrawElems();
319  SendExposeEvent();
320 }
321 
322 static void KeyMPressed()
323 {
324  vssol3d -> ToggleCPDrawMesh();
325  SendExposeEvent();
326 }
327 
328 static void KeyEPressed()
329 {
330  vssol3d -> ToggleCPDrawElems();
331  SendExposeEvent();
332 }
333 
334 static void KeyFPressed()
335 {
336  vssol3d -> ToggleShading();
337  SendExposeEvent();
338 }
339 
340 static void KeyoPressed(GLenum state)
341 {
342  if (state & KMOD_CTRL)
343  {
344  vssol3d -> ToggleDrawOrdering();
345  vssol3d -> PrepareOrderingCurve();
346  SendExposeEvent();
347  }
348  else
349  {
350  if (vssol3d -> TimesToRefine < 32)
351  {
352  cout << "Subdivision factor = " << ++vssol3d->TimesToRefine << endl;
353  if (vssol3d -> GetShading() == 2)
354  {
355  vssol3d->DoAutoscale(false);
356  vssol3d -> Prepare();
357  vssol3d -> PrepareLines();
358  vssol3d -> CPPrepare();
359  vssol3d -> PrepareLevelSurf();
360  SendExposeEvent();
361  }
362  }
363  }
364 }
365 
366 static void KeyOPressed()
367 {
368  if (vssol3d -> TimesToRefine > 1)
369  {
370  cout << "Subdivision factor = " << --vssol3d->TimesToRefine << endl;
371  if (vssol3d -> GetShading() == 2)
372  {
373  vssol3d->DoAutoscale(false);
374  vssol3d -> Prepare();
375  vssol3d -> PrepareLines();
376  vssol3d -> CPPrepare();
377  vssol3d -> PrepareLevelSurf();
378  SendExposeEvent();
379  }
380  }
381 }
382 
383 static void KeywPressed()
384 {
385  if (vssol3d -> GetShading() == 2)
386  {
387  if ( fabs(vssol3d -> FaceShiftScale += 0.01) < 0.001 )
388  {
389  vssol3d -> FaceShiftScale = 0.0;
390  }
391  cout << "New Shift Scale: " << vssol3d -> FaceShiftScale
392  << endl;
393  vssol3d -> Prepare();
394  vssol3d -> PrepareLines();
395  vssol3d -> CPPrepare();
396  SendExposeEvent();
397  }
398 }
399 
400 static void KeyWPressed()
401 {
402  if (vssol3d -> GetShading() == 2)
403  {
404  if ( fabs(vssol3d -> FaceShiftScale -= 0.01) < 0.001 )
405  {
406  vssol3d -> FaceShiftScale = 0.0;
407  }
408  cout << "New Shift Scale: " << vssol3d -> FaceShiftScale
409  << endl;
410  vssol3d -> Prepare();
411  vssol3d -> PrepareLines();
412  vssol3d -> CPPrepare();
413  SendExposeEvent();
414  }
415 }
416 
417 static void KeyuPressed()
418 {
419  vssol3d -> MoveLevelSurf(+1);
420  SendExposeEvent();
421 }
422 
423 static void KeyUPressed()
424 {
425  vssol3d -> MoveLevelSurf(-1);
426  SendExposeEvent();
427 }
428 
429 static void KeyvPressed()
430 {
431  vssol3d -> NumberOfLevelSurf(+1);
432  SendExposeEvent();
433 }
434 
435 static void KeyVPressed()
436 {
437  vssol3d -> NumberOfLevelSurf(-1);
438  SendExposeEvent();
439 }
440 
441 static int magic_key_pressed = 0;
442 void ToggleMagicKey()
443 {
444  magic_key_pressed = 1-magic_key_pressed;
445 }
446 
447 static void KeyF3Pressed(GLenum state)
448 {
449  if (state & KMOD_CTRL)
450  {
451  if (vssol3d->cut_lambda <= 0.95)
452  {
453  vssol3d->cut_lambda += 0.05;
454  vssol3d->cut_updated = false;
455  }
456  if (fabs(vssol3d->cut_lambda-1.0) < 1e-12) // snap to 1
457  {
458  vssol3d->cut_lambda = 1.0;
459  }
460  vssol3d->Prepare();
461  SendExposeEvent();
462  }
463  else
464  {
465  if (vssol3d->GetShading() == 2)
466  {
467  if (vssol3d->GetMesh()->Dimension() == 3 && vssol3d->bdrc.Width() == 0)
468  {
469  vssol3d->ComputeBdrAttrCenter();
470  }
471  if (vssol3d->GetMesh()->Dimension() == 2 && vssol3d->matc.Width() == 0)
472  {
473  vssol3d->ComputeElemAttrCenter();
474  }
475  vssol3d->shrink *= 0.9;
476  if (magic_key_pressed)
477  {
478  vssol3d -> Scale(1.11111111111111111111111);
479  }
480  SendExposeEvent();
481  vssol3d->Prepare();
482  vssol3d->PrepareLines();
483  SendExposeEvent();
484  }
485  }
486 }
487 
488 static void KeyF4Pressed(GLenum state)
489 {
490  if (state & KMOD_CTRL)
491  {
492  if (vssol3d->cut_lambda >= 0.05)
493  {
494  vssol3d->cut_lambda -= 0.05;
495  vssol3d->cut_updated = false;
496  }
497  if (fabs(vssol3d->cut_lambda-0.0) < 1e-12) // snap to 0
498  {
499  vssol3d->cut_lambda = 0.0;
500  }
501  vssol3d->Prepare();
502  SendExposeEvent();
503  }
504  else
505  {
506  if (vssol3d->GetShading() == 2)
507  {
508  if (vssol3d->GetMesh()->Dimension() == 3 && vssol3d->bdrc.Width() == 0)
509  {
510  vssol3d->ComputeBdrAttrCenter();
511  }
512  if (vssol3d->GetMesh()->Dimension() == 2 && vssol3d->matc.Width() == 0)
513  {
514  vssol3d->ComputeElemAttrCenter();
515  }
516  vssol3d->shrink *= 1.11111111111111111111111;
517  if (magic_key_pressed)
518  {
519  vssol3d -> Scale(0.9);
520  }
521  vssol3d->Prepare();
522  vssol3d->PrepareLines();
523  SendExposeEvent();
524  }
525  }
526 }
527 
528 static void KeyF11Pressed()
529 {
530  if (vssol3d->GetShading() == 2)
531  {
532  if (vssol3d->matc.Width() == 0)
533  {
534  vssol3d->ComputeElemAttrCenter();
535  }
536  vssol3d->shrinkmat *= 0.9;
537  if (magic_key_pressed)
538  {
539  vssol3d -> Scale(1.11111111111111111111111);
540  }
541  vssol3d->Prepare();
542  vssol3d->PrepareLines();
543  SendExposeEvent();
544  }
545 }
546 
547 static void KeyF12Pressed()
548 {
549  if (vssol3d->GetShading() == 2)
550  {
551  if (vssol3d->matc.Width() == 0)
552  {
553  vssol3d->ComputeElemAttrCenter();
554  }
555  vssol3d->shrinkmat *= 1.11111111111111111111111;
556  if (magic_key_pressed)
557  {
558  vssol3d -> Scale(0.9);
559  }
560  vssol3d->Prepare();
561  vssol3d->PrepareLines();
562  SendExposeEvent();
563  }
564 }
565 
566 static void KeyF8Pressed()
567 {
568  Mesh &mesh = *vssol3d->GetMesh();
569  int dim = mesh.Dimension();
570  const Array<int> &all_attr = ((dim == 3) ? mesh.bdr_attributes :
571  mesh.attributes);
572  Array<int> &attr_marker = vssol3d->bdr_attr_to_show;
573  int attr;
574  Array<int> attr_list(&attr, 1);
575 
576  cout << ((dim == 3) ? "Bdr a" : "A") << "ttributes ON: ";
577  for (int i = 0; i < all_attr.Size(); i++)
578  if (attr_marker[all_attr[i]-1])
579  {
580  cout << " " << all_attr[i];
581  }
582  cout << endl;
583 
584  cout << ((dim == 3) ? "Bdr a" : "A") << "ttribute to toggle : " << flush;
585  cin >> attr;
586  vssol3d->ToggleAttributes(attr_list);
587  SendExposeEvent();
588 }
589 
590 static void KeyF9Pressed()
591 {
592  Mesh &mesh = *vssol3d->GetMesh();
593  int dim = mesh.Dimension();
594  const Array<int> &attr_list = ((dim == 3) ? mesh.bdr_attributes :
595  mesh.attributes);
596  Array<int> &attr_marker = vssol3d->bdr_attr_to_show;
597  int attr, n, j, i;
598 
599  if (attr_list.Size() == 0)
600  {
601  return;
602  }
603  for (i = j = n = 0; i < attr_list.Size(); i++)
604  if (attr_marker[attr_list[i]-1])
605  {
606  j = i, n++;
607  }
608  if (n == 1)
609  {
610  j = (j + 1) % (attr_list.Size() + 1);
611  }
612  else
613  {
614  j = 0;
615  }
616  if (j == attr_list.Size())
617  {
618  attr_marker = 1;
619  cout << "Showing all " << ((dim == 3) ? "bdr " : "") << "attributes "
620  << endl;
621  }
622  else
623  {
624  attr = attr_list[j];
625  attr_marker = 0;
626  attr_marker[attr-1] = 1;
627  cout << "Showing " << ((dim == 3) ? "bdr " : "") << "attribute "
628  << attr << endl;
629  }
630  vssol3d -> PrepareLines();
631  vssol3d -> Prepare();
632  SendExposeEvent();
633 }
634 
635 static void KeyF10Pressed()
636 {
637  Mesh &mesh = *vssol3d->GetMesh();
638  int dim = mesh.Dimension();
639  const Array<int> &attr_list = ((dim == 3) ? mesh.bdr_attributes :
640  mesh.attributes);
641  Array<int> &attr_marker = vssol3d->bdr_attr_to_show;
642  int attr, n, j, i;
643 
644  if (attr_list.Size() == 0)
645  {
646  return;
647  }
648  n = 0;
649  for (i = j = n = 0; i < attr_list.Size(); i++)
650  if (attr_marker[attr_list[i]-1])
651  {
652  j = i, n++;
653  }
654  if (n == 1)
655  {
656  j = (j + attr_list.Size()) % (attr_list.Size() + 1);
657  }
658  else
659  {
660  j = attr_list.Size() - 1;
661  }
662  if (j == attr_list.Size())
663  {
664  attr_marker = 1;
665  cout << "Showing all " << ((dim == 3) ? "bdr " : "") << "attributes "
666  << endl;
667  }
668  else
669  {
670  attr = attr_list[j];
671  attr_marker = 0;
672  attr_marker[attr-1] = 1;
673  cout << "Showing " << ((dim == 3) ? "bdr " : "") << "attribute "
674  << attr << endl;
675  }
676  vssol3d -> PrepareLines();
677  vssol3d -> Prepare();
678  SendExposeEvent();
679 }
680 
681 VisualizationSceneSolution3d::VisualizationSceneSolution3d()
682 {}
683 
684 VisualizationSceneSolution3d::VisualizationSceneSolution3d(Mesh &m, Vector &s)
685 {
686  mesh = &m;
687  sol = &s;
688  GridF = NULL;
689 
690  Init();
691 }
692 
693 
694 void VisualizationSceneSolution3d::Init()
695 {
696  vssol3d = this;
697 
698  cplane = 0;
699  cp_drawmesh = 0; cp_drawelems = 1;
700  drawlsurf = 0;
701  cp_algo = 0;
702 
703  drawelems = shading = 1;
704  drawmesh = 0;
705  draworder = 0;
706  scaling = 0;
707 
708  shrink = 1.0;
709  shrinkmat = 1.0;
710  bdrc.SetSize(3,0);
711  matc.SetSize(3,0);
712 
713  TimesToRefine = 1;
714  FaceShiftScale = 0.0;
715 
716  if (mesh->Dimension() == 3)
717  {
718  if (!mesh->bdr_attributes.Size())
719  {
720  MFEM_WARNING("Missing boundary attributes!");
721  }
722  bdr_attr_to_show.SetSize(mesh->bdr_attributes.Size() > 0 ?
723  mesh->bdr_attributes.Max() : 0);
724  }
725  else
726  {
727  if (!mesh->attributes.Size())
728  {
729  MFEM_WARNING("Missing element attributes!");
730  }
731  bdr_attr_to_show.SetSize(mesh->attributes.Size() > 0 ?
732  mesh->attributes.Max() : 0);
733  }
734  bdr_attr_to_show = 1;
735 
736  VisualizationSceneScalarData::Init(); // calls FindNewBox
737 
738  FindNewValueRange(false);
739 
740  node_pos = new double[mesh->GetNV()];
741 
742  palette.SetIndex(12); // use the 'vivid' palette in 3D
743 
744  double eps = 1e-6; // move the cutting plane a bit to avoid artifacts
745  CuttingPlane = new Plane(-1.0,0.0,0.0,(0.5-eps)*bb.x[0]+(0.5+eps)*bb.x[1]);
746 
747  nlevels = 1;
748 
749  FindNodePos();
750 
751  // static int init = 0;
752  // if (!init)
753  {
754  // init = 1;
755 
756  wnd->setOnKeyDown('m', KeymPressed);
757  wnd->setOnKeyDown('M', KeyMPressed);
758 
759  wnd->setOnKeyDown('e', KeyePressed);
760  wnd->setOnKeyDown('E', KeyEPressed);
761 
762  wnd->setOnKeyDown('f', KeyFPressed);
763  // wnd->setOnKeyDown('F', KeyFPressed);
764 
765  wnd->setOnKeyDown('i', KeyiPressed);
766  wnd->setOnKeyDown('I', KeyIPressed);
767 
768  wnd->setOnKeyDown('o', KeyoPressed);
769  wnd->setOnKeyDown('O', KeyOPressed);
770 
771  wnd->setOnKeyDown('w', KeywPressed);
772  wnd->setOnKeyDown('W', KeyWPressed);
773 
774  wnd->setOnKeyDown('x', KeyxPressed);
775  wnd->setOnKeyDown('X', KeyXPressed);
776 
777  wnd->setOnKeyDown('y', KeyyPressed);
778  wnd->setOnKeyDown('Y', KeyYPressed);
779 
780  wnd->setOnKeyDown('z', KeyzPressed);
781  wnd->setOnKeyDown('Z', KeyZPressed);
782 
783  wnd->setOnKeyDown('u', KeyuPressed);
784  wnd->setOnKeyDown('U', KeyUPressed);
785 
786  wnd->setOnKeyDown('v', KeyvPressed);
787  wnd->setOnKeyDown('V', KeyVPressed);
788 
789  wnd->setOnKeyDown(SDLK_F3, KeyF3Pressed);
790  wnd->setOnKeyDown(SDLK_F4, KeyF4Pressed);
791  wnd->setOnKeyDown(SDLK_NUMLOCKCLEAR, ToggleMagicKey);
792 
793  wnd->setOnKeyDown(SDLK_F8, KeyF8Pressed);
794  wnd->setOnKeyDown(SDLK_F9, KeyF9Pressed);
795  wnd->setOnKeyDown(SDLK_F10, KeyF10Pressed);
796 
797  wnd->setOnKeyDown(SDLK_F11, KeyF11Pressed);
798  wnd->setOnKeyDown(SDLK_F12, KeyF12Pressed);
799  }
800  Prepare();
801  PrepareLines();
802  CPPrepare();
803  PrepareLevelSurf();
804  PrepareOrderingCurve();
805 }
806 
807 VisualizationSceneSolution3d::~VisualizationSceneSolution3d()
808 {
809  delete [] node_pos;
810 }
811 
812 void VisualizationSceneSolution3d::NewMeshAndSolution(
813  Mesh *new_m, Vector *new_sol, GridFunction *new_u)
814 {
815  if (mesh->GetNV() != new_m->GetNV())
816  {
817  delete [] node_pos;
818  node_pos = new double[new_m->GetNV()];
819  }
820  // If the number of surface elements changes, recompute the refinement factor
821  if (mesh->Dimension() != new_m->Dimension() ||
822  (mesh->Dimension() == 2 && mesh->GetNE() != new_m->GetNE()) ||
823  (mesh->Dimension() == 3 && mesh->GetNBE() != new_m->GetNBE()))
824  {
825  mesh = new_m;
826  int ref = GetAutoRefineFactor();
827  if (TimesToRefine != ref)
828  {
829  TimesToRefine = ref;
830  cout << "Subdivision factor = " << TimesToRefine << endl;
831  }
832  }
833  mesh = new_m;
834  sol = new_sol;
835  GridF = new_u;
836  FindNodePos();
837 
838  DoAutoscale(false);
839 
840  Prepare();
841  PrepareLines();
842  CPPrepare();
843  PrepareLevelSurf();
844  PrepareOrderingCurve();
845 }
846 
847 void VisualizationSceneSolution3d::SetShading(int s, bool print)
848 {
849  if (shading == s || s < 0)
850  {
851  return;
852  }
853 
854  if (s > 2 || (GridF == NULL && s > 1))
855  {
856  return;
857  }
858  int os = shading;
859  shading = s;
860  if (GridF != NULL && (s == 2 || os == 2))
861  {
862  DoAutoscale(false);
863  PrepareLines();
864  CPPrepare();
865  }
866  Prepare();
867  PrepareLevelSurf();
868 
869  static const char *shading_type[3] =
870  {"flat", "smooth", "non-conforming (with subdivision)"};
871  if (print)
872  {
873  cout << "Shading type : " << shading_type[shading] << endl;
874  }
875 }
876 
877 void VisualizationSceneSolution3d::ToggleShading()
878 {
879  if (GridF)
880  {
881  SetShading((shading+1)%3, true);
882  }
883  else
884  {
885  SetShading(1-shading, true);
886  }
887 }
888 
889 void VisualizationSceneSolution3d::SetRefineFactors(int f, int ignored)
890 {
891  if (TimesToRefine == f || f < 1)
892  {
893  return;
894  }
895 
896  TimesToRefine = f;
897 
898  if (shading == 2)
899  {
900  DoAutoscale(false);
901  Prepare();
902  PrepareLines();
903  CPPrepare();
904  PrepareOrderingCurve();
905  }
906 }
907 
908 int VisualizationSceneSolution3d::GetAutoRefineFactor()
909 {
910  int ne = mesh->GetNBE(), ref = 1;
911  if (mesh->Dimension() == 2)
912  {
913  ne = mesh->GetNE();
914  }
915 
916  while (ref < auto_ref_max && ne*(ref+1)*(ref+1) <= auto_ref_max_surf_elem)
917  {
918  ref++;
919  }
920 
921  return ref;
922 }
923 
924 void VisualizationSceneSolution3d::AutoRefine()
925 {
926  int ref = GetAutoRefineFactor();
927 
928  cout << "Subdivision factor = " << ref << endl;
929 
930  SetRefineFactors(ref, 1);
931 }
932 
933 void VisualizationSceneSolution3d::ToggleAttributes(Array<int> &attr_list)
934 {
935  int dim = mesh->Dimension();
936  Array<int> &attr_marker = bdr_attr_to_show;
937 
938  for (int i = 0; i < attr_list.Size(); i++)
939  {
940  int attr = attr_list[i];
941  if (attr < 1)
942  {
943  cout << "Hiding all" << ((dim == 3) ? " bdr" : "") << " attributes."
944  << endl;
945  attr_marker = 0;
946  }
947  else if (attr > attr_marker.Size())
948  {
949  cout << "Showing all" << ((dim == 3) ? " bdr" : "") << " attributes."
950  << endl;
951  attr_marker = 1;
952  }
953  else
954  {
955  attr_marker[attr-1] = !attr_marker[attr-1];
956  }
957  }
958  PrepareLines();
959  Prepare();
960 }
961 
962 void VisualizationSceneSolution3d::FindNewBox(bool prepare)
963 {
964  int nv = mesh -> GetNV();
965 
966  double *coord = mesh->GetVertex(0);
967 
968  bb.x[0] = bb.x[1] = coord[0];
969  bb.y[0] = bb.y[1] = coord[1];
970  bb.z[0] = bb.z[1] = coord[2];
971 
972  for (int i = 1; i < nv; i++)
973  {
974  coord = mesh->GetVertex(i);
975  if (coord[0] < bb.x[0]) { bb.x[0] = coord[0]; }
976  if (coord[1] < bb.y[0]) { bb.y[0] = coord[1]; }
977  if (coord[2] < bb.z[0]) { bb.z[0] = coord[2]; }
978  if (coord[0] > bb.x[1]) { bb.x[1] = coord[0]; }
979  if (coord[1] > bb.y[1]) { bb.y[1] = coord[1]; }
980  if (coord[2] > bb.z[1]) { bb.z[1] = coord[2]; }
981  }
982 
983  if (shading == 2)
984  {
985  int dim = mesh->Dimension();
986  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
987  int fn, fo;
988  DenseMatrix pointmat;
989  RefinedGeometry *RefG;
990  IntegrationRule eir;
991  FaceElementTransformations *Tr;
992  ElementTransformation *T;
993 
994  for (int i = 0; i < ne; i++)
995  {
996  if (dim == 3)
997  {
998  mesh->GetBdrElementFace(i, &fn, &fo);
999  RefG = GLVisGeometryRefiner.Refine(mesh->GetFaceBaseGeometry(fn),
1000  TimesToRefine);
1001  Tr = mesh->GetFaceElementTransformations(fn, 5);
1002  eir.SetSize(RefG->RefPts.GetNPoints());
1003  Tr->Loc1.Transform(RefG->RefPts, eir);
1004  Tr->Elem1->Transform(eir, pointmat);
1005  }
1006  else
1007  {
1008  T = mesh->GetElementTransformation(i);
1009  RefG = GLVisGeometryRefiner.Refine(mesh->GetElementBaseGeometry(i),
1010  TimesToRefine);
1011  T->Transform(RefG->RefPts, pointmat);
1012  }
1013  for (int j = 0; j < pointmat.Width(); j++)
1014  {
1015  if (pointmat(0,j) < bb.x[0]) { bb.x[0] = pointmat(0,j); }
1016  if (pointmat(1,j) < bb.y[0]) { bb.y[0] = pointmat(1,j); }
1017  if (pointmat(2,j) < bb.z[0]) { bb.z[0] = pointmat(2,j); }
1018  if (pointmat(0,j) > bb.x[1]) { bb.x[1] = pointmat(0,j); }
1019  if (pointmat(1,j) > bb.y[1]) { bb.y[1] = pointmat(1,j); }
1020  if (pointmat(2,j) > bb.z[1]) { bb.z[1] = pointmat(2,j); }
1021  }
1022  }
1023  }
1024 
1025  UpdateBoundingBox();
1026 }
1027 
1028 void VisualizationSceneSolution3d::FindNewValueRange(bool prepare)
1029 {
1030  if (shading < 2)
1031  {
1032  minv = sol->Min();
1033  maxv = sol->Max();
1034  }
1035  else
1036  {
1037  minv = GridF->Min();
1038  maxv = GridF->Max();
1039  }
1040  FixValueRange();
1041  UpdateValueRange(prepare);
1042 }
1043 
1044 void VisualizationSceneSolution3d::EventUpdateColors()
1045 {
1046  Prepare();
1047  PrepareCuttingPlane();
1048  PrepareLevelSurf();
1049  PrepareOrderingCurve();
1050  if (shading == 2 && drawmesh != 0 && FaceShiftScale != 0.0)
1051  {
1052  PrepareLines();
1053  }
1054 }
1055 
1056 void VisualizationSceneSolution3d::UpdateValueRange(bool prepare)
1057 {
1058  logscale = logscale && LogscaleRange();
1059  palette.SetUseLogscale(logscale);
1060  SetLogA();
1061  SetLevelLines(minv, maxv, nl);
1062  if (prepare)
1063  {
1064  UpdateLevelLines();
1065  EventUpdateColors();
1066  }
1067 }
1068 
1069 void VisualizationSceneSolution3d::FindNodePos()
1070 {
1071  int i, nnodes = mesh -> GetNV();
1072 
1073  for (i = 0; i < nnodes; i++)
1074  {
1075  node_pos[i] = CuttingPlane -> Transform (mesh -> GetVertex (i));
1076  }
1077 }
1078 
1079 void VisualizationSceneSolution3d::ToggleDrawMesh()
1080 {
1081  drawmesh = (drawmesh+1)%3;
1082  PrepareLines();
1083 }
1084 
1085 void VisualizationSceneSolution3d::ToggleCuttingPlane()
1086 {
1087  if (cplane == 2 && cp_drawmesh == 3)
1088  {
1089  cp_drawmesh = 2;
1090  }
1091 
1092  cplane = (cplane+1)%3;
1093 #ifdef GLVIS_DEBUG
1094  cout << "cplane = " << cplane << endl;
1095 #endif
1096  CPPrepare();
1097  if (cplane == 0 || cplane == 2)
1098  {
1099  Prepare();
1100  PrepareLines();
1101  PrepareOrderingCurve();
1102  }
1103 }
1104 
1105 void VisualizationSceneSolution3d::ToggleCPDrawElems()
1106 {
1107  cp_drawelems = 1-cp_drawelems;
1108  PrepareCuttingPlane();
1109 }
1110 
1111 void VisualizationSceneSolution3d::ToggleCPDrawMesh()
1112 {
1113  if (cplane == 1)
1114  {
1115  cp_drawmesh = (cp_drawmesh+1)%3;
1116  }
1117  else if (cplane == 2)
1118  {
1119  cp_drawmesh = (cp_drawmesh+1)%4;
1120  }
1121  PrepareCuttingPlaneLines();
1122 }
1123 
1124 void VisualizationSceneSolution3d::ToggleCPAlgorithm()
1125 {
1126  cp_algo = (cp_algo+1)%2;
1127  if (shading == 2 && cplane == 1)
1128  {
1129  CPPrepare();
1130  }
1131 }
1132 
1133 void VisualizationSceneSolution3d::MoveLevelSurf(int move)
1134 {
1135  drawlsurf += move;
1136  if (drawlsurf < 0)
1137  {
1138  drawlsurf = 0;
1139  }
1140  if (drawlsurf > 49)
1141  {
1142  drawlsurf = 49;
1143  }
1144  PrepareLevelSurf();
1145 }
1146 
1147 void VisualizationSceneSolution3d::NumberOfLevelSurf(int c)
1148 {
1149  nlevels += c;
1150  if (nlevels < 1)
1151  {
1152  nlevels = 1;
1153  }
1154  PrepareLevelSurf();
1155 }
1156 
1157 void VisualizationSceneSolution3d::GetFaceNormals(
1158  const int FaceNo, const int side, const IntegrationRule &ir,
1159  DenseMatrix &normals)
1160 {
1161  // match the way GridFunction::GetFaceValues works
1162  double aJInv[9], alnor[3];
1163  DenseMatrix JInv(aJInv, 3, 3);
1164  Vector lnor(alnor, 3), nr;
1165  IntegrationRule eir(ir.GetNPoints());
1166  FaceElementTransformations *Tr;
1167  normals.SetSize(3, ir.GetNPoints());
1168  ElementTransformation *ETr;
1169  IntegrationPointTransformation *LTr;
1170  if (side == 0)
1171  {
1172  Tr = mesh->GetFaceElementTransformations(FaceNo, 5);
1173  ETr = Tr->Elem1;
1174  LTr = &Tr->Loc1;
1175  }
1176  else
1177  {
1178  Tr = mesh->GetFaceElementTransformations(FaceNo, 10);
1179  ETr = Tr->Elem2;
1180  LTr = &Tr->Loc2;
1181  }
1182  LTr->Transform(ir, eir);
1183  for (int i = 0; i < normals.Width(); i++)
1184  {
1185  LTr->Transf.SetIntPoint(&ir.IntPoint(i));
1186  CalcOrtho(LTr->Transf.Jacobian(), lnor);
1187  ETr->SetIntPoint(&eir.IntPoint(i));
1188  const DenseMatrix &Jac = ETr->Jacobian();
1189  CalcInverse(Jac, JInv);
1190  normals.GetColumnReference(i, nr);
1191  JInv.MultTranspose(lnor, nr);
1192  }
1193  if (side)
1194  {
1195  normals *= -1.;
1196  }
1197  JInv.ClearExternalData();
1198 }
1199 
1200 void VisualizationSceneSolution3d::DrawRefinedSurf(
1201  int n, double *points, int elem, int func, int part)
1202 {
1203  int i, j;
1204  RefinedGeometry *RefG;
1205  IntegrationPointTransformation ip_transf;
1206 
1207  switch (n)
1208  {
1209  case 3:
1210  RefG = GLVisGeometryRefiner.Refine(Geometry::TRIANGLE, TimesToRefine);
1211  ip_transf.Transf.SetFE (&TriangleFE);
1212  break;
1213  case 4:
1214  RefG = GLVisGeometryRefiner.Refine(Geometry::SQUARE, TimesToRefine);
1215  ip_transf.Transf.SetFE (&QuadrilateralFE);
1216  break;
1217  case 5:
1218  DrawRefinedSurf (3, points, elem, func, 0); // draw (0,1,2)
1219  for (i = 0; i < 3; i++)
1220  {
1221  points[4+i] = points[i]; // move point 0 to point 1
1222  }
1223  DrawRefinedSurf (4, points+4, elem, func, 1); // draw (0,2,3,4)
1224  return;
1225  case 6:
1226  DrawRefinedSurf (4, points, elem, func, 0); // draw (0,1,2,3)
1227  for (i = 0; i < 3; i++)
1228  {
1229  points[8+i] = points[i]; // move point 0 to point 2
1230  }
1231  DrawRefinedSurf (4, points+8, elem, func, 1); // draw (0,3,4,5)
1232  return;
1233  default:
1234  return;
1235  }
1236  DenseMatrix &pm = ip_transf.Transf.GetPointMat();
1237  pm.SetSize (3, n);
1238  for (i = 0; i < n; i++)
1239  for (j = 0; j < 3; j++)
1240  {
1241  pm(j,i) = points[4*i+j];
1242  }
1243  IntegrationRule tir (RefG -> RefPts.GetNPoints());
1244  ip_transf.Transform (RefG -> RefPts, tir);
1245  DenseMatrix pointmat;
1246  Vector values;
1247  GridF -> GetValues (elem, tir, values, pointmat);
1248 
1249  LiftRefinedSurf (n, pointmat, values, NULL);
1250 
1251  switch (func)
1252  {
1253  case 1:
1254  DrawRefinedSurf (n, pointmat, values, RefG->RefGeoms);
1255  break;
1256 
1257  case 2:
1258  DrawRefinedSurfEdges (n, pointmat, values, RefG->RefEdges, part);
1259  break;
1260 
1261  case 3:
1262  DrawRefinedSurfLevelLines (n, pointmat, values, RefG->RefGeoms);
1263  break;
1264  }
1265 }
1266 
1267 void VisualizationSceneSolution3d::LiftRefinedSurf(
1268  int n, DenseMatrix &pointmat, Vector &values, int *RG)
1269 {
1270  int i, j;
1271 
1272  if (FaceShiftScale == 0.0)
1273  {
1274  return;
1275  }
1276 
1277  double norm[3];
1278 
1279  if (RG == NULL) // use the normal vector of the cutting plane
1280  {
1281  double *eqn = CuttingPlane -> Equation();
1282  norm[0] = -eqn[0]; norm[1] = -eqn[1]; norm[2] = -eqn[2];
1283  }
1284  else
1285  {
1286  // use the normal vector to the polygon defined by
1287  // the points with indexes given by the first n integers in RG
1288  double pts[4][3];
1289  for (i = 0; i < n; i++)
1290  for (j = 0; j < 3; j++)
1291  {
1292  pts[i][j] = pointmat(j,RG[i]);
1293  }
1294  if (n > 3)
1295  {
1296  j = Compute3DUnitNormal (pts[0], pts[1], pts[2], pts[3], norm);
1297  }
1298  else
1299  {
1300  j = Compute3DUnitNormal (pts[0], pts[1], pts[2], norm);
1301  }
1302  if (j)
1303  {
1304  // could not compute normal
1305  cerr << "WARNING: VisualizationSceneSolution3d::LiftRefinedSurf"
1306  << endl;
1307  return;
1308  }
1309  }
1310 
1311  double bbox_diam = sqrt ( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
1312  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
1313  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
1314  double sc = FaceShiftScale * bbox_diam;
1315 
1316  for (i = 0; i < pointmat.Width(); i++)
1317  {
1318  double val = sc * (values(i) - minv) / (maxv - minv);
1319  for (j = 0; j < 3; j++)
1320  {
1321  pointmat(j, i) += val * norm[j];
1322  }
1323  }
1324 }
1325 
1326 void VisualizationSceneSolution3d::DrawRefinedSurf(
1327  int n, DenseMatrix &pointmat, Vector &values, Array<int> &RefGeoms)
1328 {
1329  double norm[3], pts[4][3];
1330  gl3::GlBuilder draw = cplane_buf.createBuilder();
1331 
1332  for (int i = 0; i < RefGeoms.Size()/n; i++)
1333  {
1334  int *RG = &(RefGeoms[i*n]);
1335  int j;
1336  for (j = 0; j < n; j++)
1337  for (int l = 0; l < 3; l++)
1338  {
1339  pts[j][l] = pointmat(l, RG[j]);
1340  }
1341  if (n > 3)
1342  {
1343  j = Compute3DUnitNormal (pts[0], pts[1], pts[2], pts[3], norm);
1344  }
1345  else
1346  {
1347  j = Compute3DUnitNormal (pts[0], pts[1], pts[2], norm);
1348  }
1349  if (!j)
1350  {
1351  draw.glBegin (GL_POLYGON);
1352  draw.glNormal3dv (norm);
1353  for (j = 0; j < n; j++)
1354  {
1355  MySetColor (draw, values(RG[j]), minv, maxv);
1356  draw.glVertex3dv (pts[j]);
1357  }
1358  draw.glEnd();
1359  }
1360  /*
1361  else
1362  cerr << "WARNING: VisualizationSceneSolution3d::DrawRefinedSurf"
1363  << endl;
1364  */
1365  }
1366 }
1367 
1368 void VisualizationSceneSolution3d::DrawRefinedSurfEdges(
1369  int n, DenseMatrix &pointmat, Vector &values, Array<int> &RefEdges,
1370  int part)
1371 {
1372  int k, k_start, k_end;
1373 
1374  k_start = 0;
1375  k_end = RefEdges.Size();
1376  gl3::GlBuilder line = cplines_buf.createBuilder();
1377  if (part == 0)
1378  {
1379  k_end = (k_end/n) * (n-1);
1380  }
1381  if (part == 1)
1382  {
1383  k_start = (k_end/n);
1384  }
1385 
1386  line.glBegin(GL_LINES);
1387  for (k = k_start; k < k_end; k++)
1388  {
1389  int RE = RefEdges[k];
1390 
1391  line.glVertex3d (pointmat(0, RE), pointmat(1, RE),
1392  pointmat(2, RE));
1393  }
1394  line.glEnd();
1395 }
1396 
1397 void VisualizationSceneSolution3d::DrawRefinedSurfLevelLines(
1398  int n, DenseMatrix &pointmat, Vector &values, Array<int> &RefGeoms)
1399 {
1400  int j, k;
1401  int *RG;
1402  double point[4][4];
1403 
1404  gl3::GlBuilder line = cplines_buf.createBuilder();
1405 
1406  for (k = 0; k < RefGeoms.Size()/n; k++)
1407  {
1408  RG = &(RefGeoms[k*n]);
1409 
1410  for (j = 0; j < n; j++)
1411  {
1412  for (int i = 0; i < 3; i++)
1413  {
1414  point[j][i] = pointmat(i, RG[j]);
1415  }
1416  point[j][3] = values(RG[j]);
1417  }
1418  DrawPolygonLevelLines(line, point[0], n, level, false);
1419  }
1420 }
1421 
1422 void VisualizationSceneSolution3d::PrepareFlat()
1423 {
1424  int i, j;
1425  disp_buf.clear();
1426  int dim = mesh->Dimension();
1427  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
1428  DenseMatrix pointmat;
1429  Array<int> vertices;
1430  double p[4][3], c[4];
1431 
1432  for (i = 0; i < ne; i++)
1433  {
1434  if (dim == 3)
1435  {
1436  if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) { continue; }
1437 
1438  if (cplane == 2)
1439  {
1440  // for cplane == 2, get vertices of the volume element, not bdr
1441  int f, o, e1, e2;
1442  mesh->GetBdrElementFace(i, &f, &o);
1443  mesh->GetFaceElements(f, &e1, &e2);
1444  mesh->GetElementVertices(e1, vertices);
1445  }
1446  else
1447  {
1448  mesh->GetBdrElementVertices(i, vertices);
1449  }
1450  }
1451  else
1452  {
1453  if (!bdr_attr_to_show[mesh->GetAttribute(i)-1]) { continue; }
1454 
1455  mesh->GetElementVertices(i, vertices);
1456  }
1457 
1458  if (cplane == 2 && CheckPositions(vertices)) { continue; }
1459 
1460  if (dim == 3)
1461  {
1462  mesh->GetBdrPointMatrix(i, pointmat);
1463  }
1464  else
1465  {
1466  mesh->GetPointMatrix(i, pointmat);
1467  }
1468 
1469  for (j = 0; j < pointmat.Width(); j++)
1470  {
1471  p[j][0] = pointmat(0, j);
1472  p[j][1] = pointmat(1, j);
1473  p[j][2] = pointmat(2, j);
1474  c[j] = (*sol)(vertices[j]);
1475  }
1476  if (j == 3)
1477  {
1478  if (cut_lambda > 0)
1479  {
1480  DrawCutTriangle(disp_buf, p, c, minv, maxv);
1481  }
1482  else
1483  {
1484  DrawTriangle(disp_buf, p, c, minv, maxv);
1485  }
1486  }
1487  else
1488  {
1489  if (cut_lambda > 0)
1490  {
1491  DrawCutQuad(disp_buf, p, c, minv, maxv);
1492  }
1493  else
1494  {
1495  DrawQuad(disp_buf, p, c, minv, maxv);
1496  }
1497  }
1498  }
1499  updated_bufs.emplace_back(&disp_buf);
1500 }
1501 
1502 
1503 // Cut the reference square by subdividing it into 4 trapezoids with a central
1504 // square removed: (fl,fl)-(fr,fl)-(fr,fr)-(fl,fr). The input RefG corresponds
1505 // to the reference square. The value of lambda controls the cut: 0 = no cut, 1
1506 // = full cut. See keys Ctrl+F3/F4.
1507 static void CutReferenceSquare(RefinedGeometry *RefG, double lambda,
1508  IntegrationRule &RefPts, Array<int> &RefGeoms)
1509 {
1510  // lambda * vertex + (1-lambda) * center
1511  double fl = (1.0-lambda)/2.0; // left corner of the cut frame
1512  double fr = (1.0+lambda)/2.0; // right corner of the cut frame
1513 
1514  int np = RefG->RefPts.Size();
1515  RefPts.SetSize(4*np);
1516  for (int i = 0; i < np; i++)
1517  {
1518  double X = RefG->RefPts[i].x;
1519  double Y = RefG->RefPts[i].y;
1520 
1521  // First order unit square basis functions
1522  double phi3 = (1.0-X)*Y, phi2 = X*Y; // (0,1)-(1,1)
1523  double phi0 = (1.0-X)*(1.0-Y), phi1 = X*(1.0-Y); // (0,0)-(1,0)
1524 
1525  // bottom trapezoid: (0,0)-(1,1)-(fr,fl)-(fl,fl)
1526  RefPts[i].x = phi1 + fr * phi2 + fl * phi3;
1527  RefPts[i].y = fl * phi2 + fl * phi3;
1528 
1529  // right trapezoid: (fr,fl)-(1,0)-(1,1)-(fr,fr)
1530  RefPts[i+np].x = fr * phi0 + phi1 + phi2 + fr * phi3;
1531  RefPts[i+np].y = fl * phi0 + phi2 + fr * phi3;
1532 
1533  // top trapezoid: (fl,fr)-(fr,fr)-(1,1)-(0,1)
1534  RefPts[i+2*np].x = fl * phi0 + fr * phi1 + phi2;
1535  RefPts[i+2*np].y = fr * phi0 + fr * phi1 + phi2 + phi3;
1536 
1537  // left trapezoid: (0,0)-(fl,fl)-(fl,fr)-(0,1)
1538  RefPts[i+3*np].x = fl * phi1 + fl * phi2;
1539  RefPts[i+3*np].y = fl * phi1 + fr * phi2 + phi3;
1540 
1541  RefPts[i].z = RefG->RefPts[i].z;
1542  RefPts[i+np].z = RefG->RefPts[i].z;
1543  RefPts[i+2*np].z = RefG->RefPts[i].z;
1544  RefPts[i+3*np].z = RefG->RefPts[i].z;
1545  }
1546 
1547  int ne = RefG->RefGeoms.Size();
1548  RefGeoms.SetSize(4*ne);
1549  for (int i = 0; i < ne; i++)
1550  {
1551  RefGeoms[i] = RefG->RefGeoms[i];
1552  RefGeoms[i+ne] = RefG->RefGeoms[i] + np;
1553  RefGeoms[i+2*ne] = RefG->RefGeoms[i] + 2*np;
1554  RefGeoms[i+3*ne] = RefG->RefGeoms[i] + 3*np;
1555  }
1556 }
1557 
1558 // Cut the reference triangle by subdividing it into 3 trapezoids with a central
1559 // triangle removed: (fl,fl)-(fr,fl)-(fl,fr). Note that the input RefG
1560 // corresponds to a reference square, not reference triangle. The value of
1561 // lambda controls the cut: 0 = no cut, 1 = full cut. See keys Ctrl+F3/F4.
1562 static void CutReferenceTriangle(RefinedGeometry *RefG, double lambda,
1563  IntegrationRule &RefPts, Array<int> &RefGeoms)
1564 {
1565  // lambda * vertex + (1-lambda) * center
1566  double fl = (1.0-lambda)/3.0; // left corner of the cut frame
1567  double fr = (1.0+2.0*lambda)/3.0; // right corner of the cut frame
1568 
1569  int np = RefG->RefPts.Size();
1570  RefPts.SetSize(3*np);
1571  for (int i = 0; i < np; i++)
1572  {
1573  double X = RefG->RefPts[i].x;
1574  double Y = RefG->RefPts[i].y;
1575 
1576  // First order unit square basis functions
1577  double phi3 = (1.0-X)*Y, phi2 = X*Y; // (0,1)-(1,1)
1578  double phi0 = (1.0-X)*(1.0-Y), phi1 = X*(1.0-Y); // (0,0)-(1,0)
1579 
1580  // bottom trapezoid: (0,0)-(1,0)-(fr,fl)-(fl,fl)
1581  RefPts[i].x = phi1 + fr * phi2 + fl * phi3;
1582  RefPts[i].y = fl * phi2 + fl * phi3;
1583 
1584  // diagonal trapezoid: (fr,fl)-(1,0)-(0,1)-(fl,fr)
1585  RefPts[i+np].x = fr * phi0 + phi1 + fl * phi3;
1586  RefPts[i+np].y = fl * phi0 + phi2 + fr * phi3;
1587 
1588  // left trapezoid: (0,0)-(fl,fl)-(fl,fr)-(0,1)
1589  RefPts[i+2*np].x = fl * phi1 + fl * phi2;
1590  RefPts[i+2*np].y = fl * phi1 + fr * phi2 + phi3;
1591 
1592  RefPts[i].z = RefG->RefPts[i].z;
1593  RefPts[i+np].z = RefG->RefPts[i].z;
1594  RefPts[i+2*np].z = RefG->RefPts[i].z;
1595  }
1596 
1597  int ne = RefG->RefGeoms.Size();
1598  RefGeoms.SetSize(3*ne);
1599  for (int i = 0; i < ne; i++)
1600  {
1601  RefGeoms[i] = RefG->RefGeoms[i];
1602  RefGeoms[i+ne] = RefG->RefGeoms[i] + np;
1603  RefGeoms[i+2*ne] = RefG->RefGeoms[i] + 2*np;
1604  }
1605 }
1606 
1607 // Call CutReferenceTriangle and CutReferenceSquare to update the global
1608 // variables cut_TriPts, cut_TriGeoms, cut_QuadPts, cut_QuadGeoms.
1609 void CutReferenceElements(int TimesToRefine, double lambda)
1610 {
1611  RefinedGeometry *RefG =
1612  GLVisGeometryRefiner.Refine(Geometry::SQUARE, TimesToRefine);
1613  CutReferenceTriangle(RefG, lambda, cut_TriPts, cut_TriGeoms);
1614  CutReferenceSquare(RefG, lambda, cut_QuadPts, cut_QuadGeoms);
1615 }
1616 
1617 void VisualizationSceneSolution3d::PrepareFlat2()
1618 {
1619  int fn, fo, di, have_normals;
1620  double bbox_diam, vmin, vmax;
1621  disp_buf.clear();
1622 
1623  int dim = mesh->Dimension();
1624  int nbe = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
1625  DenseMatrix pointmat, normals;
1626  Vector values, normal;
1627  RefinedGeometry * RefG;
1628  Array<int> vertices;
1629  double norm[3];
1630  IsoparametricTransformation T;
1631 
1632  bbox_diam = sqrt ( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
1633  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
1634  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
1635  double sc = FaceShiftScale * bbox_diam;
1636 
1637  vmin = numeric_limits<double>::infinity();
1638  vmax = -vmin;
1639  for (int i = 0; i < nbe; i++)
1640  {
1641  int sides;
1642  switch ((dim == 3) ? mesh->GetBdrElementType(i) : mesh->GetElementType(i))
1643  {
1644  case Element::TRIANGLE:
1645  sides = 3;
1646  break;
1647 
1648  case Element::QUADRILATERAL:
1649  default:
1650  sides = 4;
1651  break;
1652  }
1653 
1654  if (dim == 3)
1655  {
1656  if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) { continue; }
1657 
1658  if (cplane == 2)
1659  {
1660  // for cplane == 2, get vertices of the volume element, not bdr
1661  int f, o, e1, e2;
1662  mesh->GetBdrElementFace(i, &f, &o);
1663  mesh->GetFaceElements(f, &e1, &e2);
1664  mesh->GetElementVertices(e1, vertices);
1665  }
1666  else
1667  {
1668  mesh->GetBdrElementVertices(i, vertices);
1669  }
1670  }
1671  else
1672  {
1673  if (!bdr_attr_to_show[mesh->GetAttribute(i)-1]) { continue; }
1674  mesh->GetElementVertices(i, vertices);
1675  }
1676 
1677  if (cplane == 2 && CheckPositions(vertices)) { continue; }
1678 
1679  if (dim == 3)
1680  {
1681  mesh -> GetBdrElementFace (i, &fn, &fo);
1682  RefG = GLVisGeometryRefiner.Refine(mesh -> GetFaceBaseGeometry (fn),
1683  TimesToRefine);
1684  if (!cut_updated)
1685  {
1686  // Update the cut version of the reference geometries
1687  CutReferenceElements(TimesToRefine, cut_lambda);
1688  cut_updated = true;
1689  }
1690 
1691  // di = GridF -> GetFaceValues (fn, 2, RefG->RefPts, values, pointmat);
1692  // this assumes the interior boundary faces are properly oriented ...
1693  di = fo % 2;
1694  if (di == 1 && !mesh->FaceIsInterior(fn))
1695  {
1696  di = 0;
1697  }
1698 
1699  IntegrationRule &RefPts = (cut_lambda > 0) ?
1700  ((sides == 3) ? cut_TriPts : cut_QuadPts) :
1701  RefG->RefPts;
1702  GridF -> GetFaceValues (fn, di, RefPts, values, pointmat);
1703  GetFaceNormals(fn, di, RefPts,normals);
1704  have_normals = 1;
1705  ShrinkPoints(pointmat, i, fn, di);
1706  }
1707  else
1708  {
1709  RefG = GLVisGeometryRefiner.Refine(mesh->GetElementBaseGeometry(i),
1710  TimesToRefine);
1711  if (!cut_updated)
1712  {
1713  // Update the cut version of the reference geometries
1714  CutReferenceElements(TimesToRefine, cut_lambda);
1715  cut_updated = true;
1716  }
1717  const IntegrationRule &ir = (cut_lambda > 0) ?
1718  ((sides == 3) ? cut_TriPts : cut_QuadPts) :
1719  RefG->RefPts;
1720  GridF->GetValues(i, ir, values, pointmat);
1721  normals.SetSize(3, values.Size());
1722  mesh->GetElementTransformation(i, &T);
1723  for (int j = 0; j < values.Size(); j++)
1724  {
1725  T.SetIntPoint(&ir.IntPoint(j));
1726  const DenseMatrix &J = T.Jacobian();
1727  normals.GetColumnReference(j, normal);
1728  CalcOrtho(J, normal);
1729  normal /= normal.Norml2();
1730  }
1731  have_normals = 1;
1732  di = 0;
1733  ShrinkPoints(pointmat, i, 0, 0);
1734  }
1735 
1736  vmin = fmin(vmin, values.Min());
1737  vmax = fmax(vmax, values.Max());
1738 
1739  // compute an average normal direction for the current face
1740  if (sc != 0.0)
1741  {
1742  for (int j = 0; j < 3; j++)
1743  {
1744  norm[j] = 0.0;
1745  }
1746  Normalize(normals);
1747  for (int k = 0; k < normals.Width(); k++)
1748  for (int j = 0; j < 3; j++)
1749  {
1750  norm[j] += normals(j, k);
1751  }
1752  Normalize(norm);
1753  for (int k = 0; k < pointmat.Width(); k++)
1754  {
1755  double val = sc * (values(k) - minv) / (maxv - minv);
1756  for (int j = 0; j < 3; j++)
1757  {
1758  pointmat(j, k) += val * norm[j];
1759  }
1760  }
1761  have_normals = 0;
1762  }
1763 
1764  have_normals = have_normals ? 2 : 0;
1765  if (di)
1766  {
1767  have_normals = -1 - have_normals;
1768  }
1769  // Comment the above lines and use the below version in order to remove
1770  // the 3D dark artifacts (indicating wrong boundary element orientation)
1771  // have_normals = have_normals ? 1 : 0;
1772 
1773  Array<int> &RefGeoms = (cut_lambda > 0) ?
1774  ((sides == 3) ? cut_TriGeoms : cut_QuadGeoms) :
1775  RefG->RefGeoms;
1776  int psides = (cut_lambda > 0) ? 4 : sides;
1777  DrawPatch(disp_buf, pointmat, values, normals, psides, RefGeoms,
1778  minv, maxv, have_normals);
1779  }
1780  updated_bufs.emplace_back(&disp_buf);
1781 
1782  cout << "VisualizationSceneSolution3d::PrepareFlat2() : [min,max] = ["
1783  << vmin << "," << vmax << "]" << endl;
1784 }
1785 
1786 void VisualizationSceneSolution3d::Prepare()
1787 {
1788  int j;
1789 
1790  if (!drawelems)
1791  {
1792  return;
1793  }
1794 
1795  switch (shading)
1796  {
1797  case 0:
1798  PrepareFlat();
1799  return;
1800  case 2:
1801  PrepareFlat2();
1802  return;
1803  default:
1804  break;
1805  }
1806 
1807  disp_buf.clear();
1808  gl3::GlBuilder poly = disp_buf.createBuilder();
1809 
1810  int dim = mesh->Dimension();
1811  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
1812  int nv = mesh -> GetNV();
1813  DenseMatrix pointmat;
1814  Array<int> vertices;
1815  double nor[3];
1816 
1817  Vector nx(nv);
1818  Vector ny(nv);
1819  Vector nz(nv);
1820 
1821  Table ba_to_be; // boundary_attribute--to--boundary_element
1822  {
1823  Table be_to_ba;
1824  be_to_ba.MakeI(ne);
1825  for (int i = 0; i < ne; i++)
1826  {
1827  be_to_ba.AddAColumnInRow(i);
1828  }
1829  be_to_ba.MakeJ();
1830  if (dim == 3)
1831  {
1832  for (int i = 0; i < ne; i++)
1833  {
1834  be_to_ba.AddConnection(i, mesh->GetBdrAttribute(i)-1);
1835  }
1836  }
1837  else
1838  {
1839  for (int i = 0; i < ne; i++)
1840  {
1841  be_to_ba.AddConnection(i, mesh->GetAttribute(i)-1);
1842  }
1843  }
1844  be_to_ba.ShiftUpI();
1845 
1846  Transpose(be_to_ba, ba_to_be);
1847  }
1848 
1849  const Array<int> &attributes =
1850  ((dim == 3) ? mesh->bdr_attributes : mesh->attributes);
1851  for (int d = 0; d < attributes.Size(); d++)
1852  {
1853  const int attr = attributes[d]-1;
1854 
1855  if (!bdr_attr_to_show[attr]) { continue; }
1856 
1857  const int nelem = ba_to_be.RowSize(attr);
1858  const int *elem = ba_to_be.GetRow(attr);
1859 
1860  for (int i = 0; i < nelem; i++)
1861  {
1862  if (dim == 3)
1863  {
1864  mesh->GetBdrElementVertices(elem[i], vertices);
1865  }
1866  else
1867  {
1868  mesh->GetElementVertices(elem[i], vertices);
1869  }
1870  for (j = 0; j < vertices.Size(); j++)
1871  {
1872  nx(vertices[j]) = ny(vertices[j]) = nz(vertices[j]) = 0.;
1873  }
1874  }
1875  for (int i = 0; i < nelem; i++)
1876  {
1877  if (dim == 3)
1878  {
1879  mesh->GetBdrPointMatrix(elem[i], pointmat);
1880  mesh->GetBdrElementVertices(elem[i], vertices);
1881  }
1882  else
1883  {
1884  mesh->GetPointMatrix(elem[i], pointmat);
1885  mesh->GetElementVertices(elem[i], vertices);
1886  }
1887 
1888  if (pointmat.Width() == 3)
1889  j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1),
1890  &pointmat(0,2), nor);
1891  else
1892  j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1),
1893  &pointmat(0,2), &pointmat(0,3), nor);
1894  if (j == 0)
1895  for (j = 0; j < pointmat.Size(); j++)
1896  {
1897  nx(vertices[j]) += nor[0];
1898  ny(vertices[j]) += nor[1];
1899  nz(vertices[j]) += nor[2];
1900  }
1901  }
1902 
1903  for (int i = 0; i < nelem; i++)
1904  {
1905  if (dim == 3)
1906  {
1907  if (cplane == 2)
1908  {
1909  // for cplane == 2, get vertices of the volume element, not bdr
1910  int f, o, e1, e2;
1911  mesh->GetBdrElementFace(i, &f, &o);
1912  mesh->GetFaceElements(f, &e1, &e2);
1913  mesh->GetElementVertices(e1, vertices);
1914 
1915  if (CheckPositions(vertices)) { continue; }
1916  }
1917  else
1918  {
1919  mesh->GetBdrElementVertices(elem[i], vertices);
1920  }
1921  }
1922  else
1923  {
1924  mesh->GetElementVertices(elem[i], vertices);
1925  }
1926 
1927  GLenum elemType = GL_NONE;
1928  switch ((dim == 3) ? mesh->GetBdrElementType(elem[i]) :
1929  mesh->GetElementType(elem[i]))
1930  {
1931  case Element::TRIANGLE:
1932  elemType = GL_TRIANGLES;
1933  break;
1934  case Element::QUADRILATERAL:
1935  elemType = GL_QUADS;
1936  break;
1937  default:
1938  MFEM_ABORT("Invalid boundary element type");
1939  break;
1940  }
1941  poly.glBegin(elemType);
1942  if (dim == 3)
1943  {
1944  mesh->GetBdrPointMatrix(elem[i], pointmat);
1945  }
1946  else
1947  {
1948  mesh->GetPointMatrix(elem[i], pointmat);
1949  }
1950 
1951  for (j = 0; j < pointmat.Size(); j++)
1952  {
1953  MySetColor(poly, (*sol)(vertices[j]), minv, maxv);
1954  poly.glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j]));
1955  poly.glVertex3dv(&pointmat(0, j));
1956  }
1957  poly.glEnd();
1958  }
1959  }
1960  updated_bufs.emplace_back(&disp_buf);
1961 }
1962 
1963 void VisualizationSceneSolution3d::PrepareLines()
1964 {
1965  if (!drawmesh)
1966  {
1967  return;
1968  }
1969 
1970  if (shading == 2)
1971  {
1972  PrepareLines2();
1973  return;
1974  }
1975 
1976  int dim = mesh->Dimension();
1977  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
1978  int i, j, k;
1979  DenseMatrix pointmat;
1980 
1981  line_buf.clear();
1982 
1983  Array<int> vertices;
1984 
1985  for (i = 0; i < ne; i++)
1986  {
1987  if (dim == 3)
1988  {
1989  if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) { continue; }
1990 
1991  if (cplane == 2)
1992  {
1993  // for cplane == 2, get vertices of the volume element, not bdr
1994  int f, o, e1, e2;
1995  mesh->GetBdrElementFace(i, &f, &o);
1996  mesh->GetFaceElements(f, &e1, &e2);
1997  mesh->GetElementVertices(e1, vertices);
1998  }
1999  else
2000  {
2001  mesh->GetBdrElementVertices(i, vertices);
2002  }
2003  }
2004  else
2005  {
2006  if (!bdr_attr_to_show[mesh->GetAttribute(i)-1]) { continue; }
2007 
2008  mesh->GetElementVertices(i, vertices);
2009  }
2010 
2011  if (cplane == 2 && CheckPositions(vertices)) { continue; }
2012 
2013  double point[4][4];
2014  if (dim == 3)
2015  {
2016  mesh->GetBdrPointMatrix(i, pointmat);
2017  }
2018  else
2019  {
2020  mesh->GetPointMatrix(i, pointmat);
2021  }
2022 
2023  gl3::GlBuilder line = line_buf.createBuilder();
2024  switch (drawmesh)
2025  {
2026  case 1:
2027  line.glBegin(GL_LINE_LOOP);
2028 
2029  for (j = 0; j < pointmat.Size(); j++)
2030  {
2031  line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j));
2032  }
2033  line.glEnd();
2034  break;
2035 
2036  case 2:
2037  for (j = 0; j < pointmat.Size(); j++)
2038  {
2039  for (k = 0; k < 3; k++)
2040  {
2041  point[j][k] = pointmat(k,j);
2042  }
2043  point[j][3] = (*sol)(vertices[j]);
2044  }
2045  DrawPolygonLevelLines(line, point[0], pointmat.Size(), level, false);
2046  break;
2047  }
2048  }
2049  updated_bufs.emplace_back(&line_buf);
2050 }
2051 
2052 void VisualizationSceneSolution3d::PrepareLines2()
2053 {
2054  int fn, fo, di = 0;
2055  double bbox_diam;
2056 
2057  line_buf.clear();
2058 
2059  int dim = mesh->Dimension();
2060  int nbe = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
2061  DenseMatrix pointmat, normals;
2062  Vector values, normal;
2063  RefinedGeometry * RefG;
2064  Array<int> vertices;
2065  IsoparametricTransformation T;
2066 
2067  bbox_diam = sqrt ( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
2068  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
2069  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
2070  double sc = FaceShiftScale * bbox_diam;
2071 
2072  for (int i = 0; i < nbe; i++)
2073  {
2074  if (dim == 3)
2075  {
2076  if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) { continue; }
2077 
2078  if (cplane == 2)
2079  {
2080  // for cplane == 2, get vertices of the volume element, not bdr
2081  int f, o, e1, e2;
2082  mesh->GetBdrElementFace(i, &f, &o);
2083  mesh->GetFaceElements(f, &e1, &e2);
2084  mesh->GetElementVertices(e1, vertices);
2085  }
2086  else
2087  {
2088  mesh->GetBdrElementVertices(i, vertices);
2089  }
2090  }
2091  else
2092  {
2093  if (!bdr_attr_to_show[mesh->GetAttribute(i)-1]) { continue; }
2094 
2095  mesh->GetElementVertices(i, vertices);
2096  }
2097 
2098  if (cplane == 2 && CheckPositions(vertices)) { continue; }
2099 
2100  if (dim == 3)
2101  {
2102  mesh -> GetBdrElementFace (i, &fn, &fo);
2103  RefG = GLVisGeometryRefiner.Refine(mesh -> GetFaceBaseGeometry (fn),
2104  TimesToRefine);
2105  // di = GridF -> GetFaceValues (fn, 2, RefG->RefPts, values, pointmat);
2106  di = fo % 2;
2107  if (di == 1 && !mesh->FaceIsInterior(fn))
2108  {
2109  di = 0;
2110  }
2111  GridF -> GetFaceValues (fn, di, RefG->RefPts, values, pointmat);
2112  ShrinkPoints(pointmat, i, fn, di);
2113  }
2114  else
2115  {
2116  RefG = GLVisGeometryRefiner.Refine(mesh->GetElementBaseGeometry(i),
2117  TimesToRefine);
2118  GridF->GetValues(i, RefG->RefPts, values, pointmat);
2119  ShrinkPoints(pointmat, i, 0, 0);
2120  }
2121 
2122  if (sc != 0.0)
2123  {
2124  if (dim == 3)
2125  {
2126  GetFaceNormals(fn, di, RefG->RefPts, normals);
2127  }
2128  else
2129  {
2130  normals.SetSize(3, values.Size());
2131  mesh->GetElementTransformation(i, &T);
2132  for (int j = 0; j < values.Size(); j++)
2133  {
2134  T.SetIntPoint(&RefG->RefPts.IntPoint(j));
2135  const DenseMatrix &J = T.Jacobian();
2136  normals.GetColumnReference(j, normal);
2137  CalcOrtho(J, normal);
2138  normal /= normal.Norml2();
2139  }
2140  }
2141  double norm[3];
2142  for (int j = 0; j < 3; j++)
2143  {
2144  norm[j] = 0.0;
2145  }
2146  for (int k = 0; k < normals.Width(); k++)
2147  {
2148  for (int j = 0; j < 3; j++)
2149  {
2150  norm[j] += normals(j, k);
2151  }
2152  }
2153  double len = sqrt(InnerProd(norm, norm));
2154  if (len > 0.0)
2155  {
2156  len = 1.0 / len;
2157  }
2158  for (int j = 0; j < 3; j++)
2159  {
2160  norm[j] *= len;
2161  }
2162  for (int k = 0; k < pointmat.Width(); k++)
2163  {
2164  double val = sc * (values(k) - minv) / (maxv - minv);
2165  for (int j = 0; j < 3; j++)
2166  {
2167  pointmat(j, k) += val * norm[j];
2168  }
2169  }
2170  }
2171  gl3::GlBuilder line = line_buf.createBuilder();
2172  if (drawmesh == 1)
2173  {
2174  Array<int> &REdges = RefG->RefEdges;
2175 
2176  line.glBegin(GL_LINES);
2177  for (int k = 0; k < REdges.Size(); k++)
2178  {
2179  line.glVertex3dv(&pointmat(0, REdges[k]));
2180  }
2181  line.glEnd();
2182  }
2183  else if (drawmesh == 2)
2184  {
2185  double point[4][4];
2186  int sides;
2187  switch ((dim == 3) ? mesh->GetBdrElementType(i) :
2188  mesh->GetElementType(i))
2189  {
2190  case Element::TRIANGLE:
2191  sides = 3;
2192  break;
2193 
2194  case Element::QUADRILATERAL:
2195  default:
2196  sides = 4;
2197  break;
2198  }
2199  for (int k = 0; k < RefG->RefGeoms.Size()/sides; k++)
2200  {
2201  int *RG = &(RefG->RefGeoms[k*sides]);
2202 
2203  for (int j = 0; j < sides; j++)
2204  {
2205  for (int ii = 0; ii < 3; ii++)
2206  {
2207  point[j][ii] = pointmat(ii, RG[j]);
2208  }
2209  point[j][3] = values(RG[j]);
2210  }
2211  DrawPolygonLevelLines(line, point[0], sides, level, false);
2212  }
2213  }
2214  }
2215  updated_bufs.emplace_back(&line_buf);
2216 }
2217 
2218 static void CutElement(const Geometry::Type geom, const int *vert_flags,
2219  const int **edge_vert_ptr, int *cut_edges,
2220  int *num_cut_edges, int *n2_cut_edges)
2221 {
2222  static const int tet_edges[12]= {0,3, 0,2, 0,1, 1,2, 1,3, 2,3};
2223  static const int pyr_edges[16] =
2224  {0,1, 1,2, 3,2, 0,3, 0,4, 1,4, 2,4, 3,4};
2225  static const int pyr_cutting[16][3] =
2226  {
2227  { 3, 4, 6}, { 9, -1, 10}, { 4, 6, 1}, {11, -1, 12},
2228  {13, -1, 14}, { 6, 1, 3}, {15, -1, 8}, { 1, 3, 4},
2229  {-1, 10, 0}, { 7, 15, -1}, {-1, 12, 2}, { 0, 9, -1},
2230  {-1, 14, 5}, { 2, 11, -1}, {-1, 8, 7}, { 5, 13, -1}
2231  };
2232  static const int pri_edges[18] =
2233  {0,1, 1,2, 2,0, 3,4, 4,5, 5,3, 0,3, 1,4, 2,5};
2234  static const int pri_cutting[18][3] =
2235  {
2236  { 3, -1, 5}, {13, 7, 14}, { 5, -1, 1}, {15, 9, 16},
2237  { 1, -1, 3}, {17, 11, 12},
2238  {14, 0, 13}, {10, -1, 8}, {16, 2, 15}, { 6, -1, 10},
2239  {12, 4, 17}, { 8, -1, 6},
2240  { 7, 14, 0}, { 4, 17, 11}, { 9, 16, 2}, { 0, 13, 7},
2241  {11, 12, 4}, { 2, 15, 9}
2242  };
2243  static const int hex_edges[24] =
2244  { 0,1, 1,2, 3,2, 0,3, 4,5, 5,6, 7,6, 4,7, 0,4, 1,5, 2,6, 3,7 };
2245  static const int hex_cutting[24][3] =
2246  {
2247  { 3, 4, 6}, {17, 9, 18}, { 4, 6, 1}, {19, 11, 20},
2248  {21, 12, 22}, { 6, 1, 3}, {23, 14, 16}, { 1, 3, 4},
2249  {18, 0, 17}, {15, 13, 10}, {20, 2, 19}, { 8, 15, 13},
2250  {10, 8, 15}, {22, 5, 21}, {13, 10, 8}, {16, 7, 23},
2251  { 9, 18, 0}, { 7, 23, 14}, {11, 20, 2}, { 0, 17, 9},
2252  {12, 22, 5}, { 2, 19, 11}, {14, 16, 7}, { 5, 21, 12}
2253  };
2254  const int *ev;
2255  int n, n2;
2256 
2257  n = n2 = 0;
2258  switch (geom)
2259  {
2260  case Geometry::TETRAHEDRON:
2261  {
2262  ev = tet_edges;
2263  for (int j = 0; j < 6; j++, ev += 2)
2264  if (vert_flags[ev[0]] != vert_flags[ev[1]])
2265  {
2266  cut_edges[n++] = j;
2267  }
2268  ev = tet_edges;
2269  }
2270  break;
2271 
2272  case Geometry::PYRAMID:
2273  {
2274  int emark[8];
2275  ev = pyr_edges;
2276  for (int j = 0; j < 8; j++, ev += 2)
2277  {
2278  emark[j] = vert_flags[ev[1]] - vert_flags[ev[0]];
2279  }
2280  do
2281  {
2282  int j;
2283  for (j = 0; j < 8; j++)
2284  {
2285  if (emark[j]) { break; }
2286  }
2287  if (j == 8)
2288  {
2289  break;
2290  }
2291  int k = 2 * j;
2292  if (emark[j] > 0)
2293  {
2294  k++;
2295  }
2296  do
2297  {
2298  int m;
2299  for (j = 0; j < 3; j++)
2300  {
2301  m = pyr_cutting[k][j];
2302  if (m >= 0)
2303  {
2304  ev = pyr_edges + 2 * (m / 2);
2305  if ((m%2 == 0 && vert_flags[ev[0]] > vert_flags[ev[1]]) ||
2306  (m%2 == 1 && vert_flags[ev[1]] > vert_flags[ev[0]]))
2307  {
2308  break;
2309  }
2310  }
2311  }
2312  cut_edges[n2++] = k/2;
2313  emark[k/2] = 0;
2314  k = m;
2315  }
2316  while (k/2 != cut_edges[n]);
2317  if (n == 0)
2318  {
2319  n = n2;
2320  }
2321  else
2322  {
2323  break;
2324  }
2325  }
2326  while (1);
2327  n2 -= n;
2328  ev = pyr_edges;
2329  }
2330  break;
2331 
2332  case Geometry::PRISM:
2333  {
2334  int emark[9];
2335  ev = pri_edges;
2336  for (int j = 0; j < 9; j++, ev += 2)
2337  {
2338  emark[j] = vert_flags[ev[1]] - vert_flags[ev[0]];
2339  }
2340  do
2341  {
2342  int j;
2343  for (j = 0; j < 9; j++)
2344  {
2345  if (emark[j]) { break; }
2346  }
2347  if (j == 9)
2348  {
2349  break;
2350  }
2351  int k = 2 * j;
2352  if (emark[j] > 0)
2353  {
2354  k++;
2355  }
2356  do
2357  {
2358  int m;
2359  for (j = 0; j < 3; j++)
2360  {
2361  m = pri_cutting[k][j];
2362  if (m >= 0)
2363  {
2364  ev = pri_edges + 2 * (m / 2);
2365  if ((m%2 == 0 && vert_flags[ev[0]] > vert_flags[ev[1]]) ||
2366  (m%2 == 1 && vert_flags[ev[1]] > vert_flags[ev[0]]))
2367  {
2368  break;
2369  }
2370  }
2371  }
2372  cut_edges[n2++] = k/2;
2373  emark[k/2] = 0;
2374  k = m;
2375  }
2376  while (k/2 != cut_edges[n]);
2377  if (n == 0)
2378  {
2379  n = n2;
2380  }
2381  else
2382  {
2383  break;
2384  }
2385  }
2386  while (1);
2387  n2 -= n;
2388  ev = pri_edges;
2389  }
2390  break;
2391 
2392  case Geometry::CUBE:
2393  {
2394  int emark[12];
2395  ev = hex_edges;
2396  for (int j = 0; j < 12; j++, ev += 2)
2397  {
2398  emark[j] = vert_flags[ev[1]] - vert_flags[ev[0]];
2399  }
2400  do
2401  {
2402  int j;
2403  for (j = 0; j < 12; j++)
2404  {
2405  if (emark[j]) { break; }
2406  }
2407  if (j == 12)
2408  {
2409  break;
2410  }
2411  int k = 2 * j;
2412  if (emark[j] > 0)
2413  {
2414  k++;
2415  }
2416  do
2417  {
2418  int m;
2419  for (j = 0; j < 3; j++)
2420  {
2421  m = hex_cutting[k][j];
2422  ev = hex_edges + 2 * (m / 2);
2423  if ((m%2 == 0 && vert_flags[ev[0]] > vert_flags[ev[1]]) ||
2424  (m%2 == 1 && vert_flags[ev[1]] > vert_flags[ev[0]]))
2425  {
2426  break;
2427  }
2428  }
2429  cut_edges[n2++] = k/2;
2430  emark[k/2] = 0;
2431  k = m;
2432  }
2433  while (k/2 != cut_edges[n]);
2434  if (n == 0)
2435  {
2436  n = n2;
2437  }
2438  else
2439  {
2440  break;
2441  }
2442  }
2443  while (1);
2444  n2 -= n;
2445  ev = hex_edges;
2446  }
2447  break;
2448 
2449  default:
2450  ev = NULL; // suppress a warning
2451  break;
2452  }
2453 
2454  *edge_vert_ptr = ev;
2455  *num_cut_edges = n;
2456  *n2_cut_edges = n2;
2457 }
2458 
2459 void VisualizationSceneSolution3d::CuttingPlaneFunc(int func)
2460 {
2461  int m, n, n2;
2462  int flag[8], cut_edges[6];
2463  const int *ev;
2464  double t, point[6][4], norm[3];
2465 
2466  DenseMatrix pointmat;
2467 
2468  Array<int> nodes;
2469  for (int i = 0; i < mesh -> GetNE(); i++)
2470  {
2471  n = n2 = 0; // n will be the number of intersection points
2472  mesh -> GetElementVertices(i, nodes);
2473  for (int j = 0; j < nodes.Size(); j++)
2474  {
2475  if (node_pos[nodes[j]] >= 0.0)
2476  {
2477  flag[j] = 1;
2478  }
2479  else
2480  {
2481  flag[j] = -1;
2482  }
2483  }
2484 
2485  CutElement(mesh->GetElementBaseGeometry(i), flag,
2486  &ev, cut_edges, &n, &n2);
2487 
2488  while (n > 2)
2489  {
2490  if (shading != 2)
2491  {
2492  mesh -> GetPointMatrix (i, pointmat);
2493  }
2494  else
2495  {
2496  const IntegrationRule *ir;
2497  ir = Geometries.GetVertices (mesh -> GetElementBaseGeometry(i));
2498  pointmat.SetSize (3, ir -> GetNPoints());
2499  for (int j = 0; j < ir -> GetNPoints(); j++)
2500  {
2501  const IntegrationPoint &ip = ir -> IntPoint (j);
2502  pointmat(0,j) = ip.x;
2503  pointmat(1,j) = ip.y;
2504  pointmat(2,j) = ip.z;
2505  }
2506  }
2507  for (int j = 0; j < n; j++)
2508  {
2509  const int *en = ev + 2*cut_edges[j];
2510  t = node_pos[ nodes[en[1]] ];
2511  t = t / ( t - node_pos[ nodes[en[0]] ] );
2512  for (int k = 0; k < 3; k++)
2513  {
2514  point[j][k] = t*pointmat(k,en[0]) + (1-t)*pointmat(k,en[1]);
2515  }
2516  point[j][3] = t*(*sol)(nodes[en[0]]) + (1-t)*(*sol)(nodes[en[1]]);
2517  }
2518 
2519  switch (func)
2520  {
2521  case 1: // PrepareCuttingPlane()
2522  {
2523  if (shading == 2)
2524  {
2525  // changes point for n > 4
2526  DrawRefinedSurf(n, point[0], i, 1);
2527  }
2528  else
2529  {
2530  m = n;
2531  int no_norm;
2532  while (1)
2533  {
2534  if (m > 3)
2535  {
2536  no_norm = Compute3DUnitNormal(point[0], point[1], point[2],
2537  point[3], norm);
2538  if (no_norm && m > 4)
2539  {
2540  for (int j = 3; j < m; j++)
2541  for (int k = 0; k < 4; k++)
2542  {
2543  point[j-2][k] = point[j][k];
2544  }
2545  m -= 2;
2546  continue;
2547  }
2548  }
2549  else
2550  no_norm = Compute3DUnitNormal(point[0], point[1], point[2],
2551  norm);
2552  break;
2553  }
2554 
2555  gl3::GlBuilder draw = cplane_buf.createBuilder();
2556  if (!no_norm)
2557  {
2558  draw.glBegin(GL_POLYGON);
2559  draw.glNormal3dv(norm);
2560  for (int j = 0; j < m; j++)
2561  {
2562  MySetColor(draw, point[j][3], minv, maxv);
2563  draw.glVertex3dv(point[j]);
2564  }
2565  draw.glEnd();
2566  }
2567  }
2568  }
2569  break;
2570 
2571  case 2: // PrepareCuttingPlaneLines() with mesh
2572  {
2573  if (shading == 2)
2574  {
2575  // changes point for n > 4
2576  DrawRefinedSurf(n, point[0], i, 2);
2577  }
2578  else
2579  {
2580  // glBegin (GL_POLYGON);
2581  gl3::GlBuilder line = cplines_buf.createBuilder();
2582  line.glBegin(GL_LINE_LOOP);
2583  for (int j = 0; j < n; j++)
2584  {
2585  line.glVertex3dv(point[j]);
2586  }
2587  line.glEnd();
2588  }
2589  }
2590  break;
2591 
2592  case 3: // PrepareCuttingPlaneLines() with level lines
2593  {
2594  if (shading == 2)
2595  {
2596  // changes point for n > 4
2597  DrawRefinedSurf(n, point[0], i, 3);
2598  }
2599  else
2600  {
2601  gl3::GlBuilder line = cplines_buf.createBuilder();
2602  DrawPolygonLevelLines(line, point[0], n, level, false);
2603  }
2604  }
2605  break;
2606  }
2607 
2608  for (int j = 0; j < n2; j++)
2609  {
2610  cut_edges[j] = cut_edges[j+n];
2611  }
2612  n = n2;
2613  n2 = 0;
2614  }
2615  }
2616 }
2617 
2618 void VisualizationSceneSolution3d::CutRefinedElement(
2619  gl3::GlDrawable& target,
2620  const DenseMatrix &verts, const Vector &vert_dist, const Vector &vals,
2621  const Geometry::Type geom, const int *elems, int num_elems, int func)
2622 {
2623  double sc = 0.0;
2624  if (FaceShiftScale != 0.0)
2625  {
2626  double bbox_diam = sqrt ( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
2627  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
2628  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
2629  sc = FaceShiftScale * bbox_diam;
2630  }
2631  const int nv = Geometry::NumVerts[geom];
2632 
2633  gl3::GlBuilder bld = target.createBuilder();
2634 
2635  for (int i = 0; i < num_elems; i++)
2636  {
2637  Geometry::Type egeom = geom;
2638  int vert_flag[8], cut_edges[8];
2639  const int *elem = elems + i*nv;
2640  const int *edge_vert;
2641  int n = 0, n2, nev = 0;
2642  for (int j = 0; j < nv; j++)
2643  {
2644  if (elem[j] < 0)
2645  {
2646  // This appears to be a tetrahedral subelement of a refined pyramid
2647  egeom = Geometry::TETRAHEDRON;
2648  break;
2649  }
2650  nev++;
2651  vert_flag[j] = (vert_dist(elem[j]) >= 0.0) ? n++, 1 : 0;
2652  }
2653  if (n == 0 || n == nev) { continue; }
2654 
2655  CutElement(egeom, vert_flag, &edge_vert, cut_edges, &n, &n2);
2656  // n = number of intersected edges
2657  // n2 = number of intersected edges, second polygon
2658 
2659  while (n > 2)
2660  {
2661  // 'pts' describe the intersecting polygon: triangle, quad, etc
2662  double pts[6][4]; // up to 6 points x (3 coordinates + 1 value)
2663 
2664  for (int j = 0; j < n; j++)
2665  {
2666  const int *ev = edge_vert + 2*cut_edges[j];
2667  double t = vert_dist(elem[ev[0]]);
2668  t = t / (t - vert_dist(elem[ev[1]]));
2669  for (int d = 0; d < 3; d++)
2670  {
2671  pts[j][d] = (1-t)*verts(d,elem[ev[0]]) + t*verts(d,elem[ev[1]]);
2672  }
2673  pts[j][3] = (1-t)*vals(elem[ev[0]]) + t*vals(elem[ev[1]]);
2674  }
2675  if (sc != 0.0)
2676  {
2677  const double *dir = CuttingPlane->Equation();
2678  for (int j = 0; j < n; j++)
2679  {
2680  // dir points into the visible side, so we add a minus to val:
2681  const double val = -sc * (pts[j][3] - minv) / (maxv - minv);
2682  for (int d = 0; d < 3; d++)
2683  {
2684  pts[j][d] += val*dir[d];
2685  }
2686  }
2687  }
2688 
2689  if (func == 0) // draw surface
2690  {
2691  double norm[3];
2692  int no_norm, m = n;
2693  while (1)
2694  {
2695  if (m > 3)
2696  {
2697  no_norm = Compute3DUnitNormal(pts[0], pts[1], pts[2], pts[3],
2698  norm);
2699  if (no_norm && m > 4)
2700  {
2701  for (int j = 3; j < m; j++)
2702  {
2703  for (int k = 0; k < 4; k++)
2704  {
2705  pts[j-2][k] = pts[j][k];
2706  }
2707  }
2708  m -= 2;
2709  continue;
2710  }
2711  }
2712  else
2713  {
2714  no_norm = Compute3DUnitNormal(pts[0], pts[1], pts[2], norm);
2715  }
2716  break;
2717  }
2718  if (!no_norm)
2719  {
2720  bld.glBegin(GL_POLYGON);
2721  bld.glNormal3dv(norm);
2722  for (int j = 0; j < m; j++)
2723  {
2724  MySetColor(bld, pts[j][3], minv, maxv);
2725  bld.glVertex3dv(pts[j]);
2726  }
2727  bld.glEnd();
2728  }
2729  }
2730  else // draw level lines
2731  {
2732  DrawPolygonLevelLines(bld, pts[0], n, level, false);
2733  }
2734 
2735  for (int j = 0; j < n2; j++)
2736  {
2737  cut_edges[j] = cut_edges[j+n];
2738  }
2739  n = n2;
2740  n2 = 0;
2741  }
2742  }
2743 }
2744 
2745 void VisualizationSceneSolution3d::CutRefinedFace(
2746  gl3::GlDrawable& target,
2747  const DenseMatrix &verts, const Vector &vert_dist, const Vector &vals,
2748  const Geometry::Type geom, const int *faces, int num_faces)
2749 {
2750  double sc = 0.0;
2751  if (FaceShiftScale != 0.0)
2752  {
2753  double bbox_diam = sqrt ( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
2754  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
2755  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
2756  sc = FaceShiftScale * bbox_diam;
2757  }
2758  const int nv = Geometry::NumVerts[geom];
2759 
2760  for (int i = 0; i < num_faces; i++)
2761  {
2762  int vert_flag[4], cut_edges[4];
2763  const int *face = faces + i*nv;
2764  int n = 0;
2765  for (int j = 0; j < nv; j++)
2766  {
2767  vert_flag[j] = (vert_dist(face[j]) >= 0.0) ? n++, 1 : 0;
2768  }
2769  if (n == 0 || n == nv) { continue; }
2770  n = 0;
2771  for (int j = 0; j < nv; j++)
2772  {
2773  const int j1 = (j+1)%nv;
2774  if (vert_flag[j] != vert_flag[j1])
2775  {
2776  cut_edges[n++] = j;
2777  }
2778  }
2779  // n = number of intersected edges (2, or 4)
2780  double pts[4][4]; // up to 4 points x (3 coordinates + 1 value)
2781  for (int j = 0; j < n; j++)
2782  {
2783  const int v0 = cut_edges[j];
2784  const int v1 = (v0+1)%nv;
2785  double t = vert_dist(face[v0]);
2786  t = t / (t - vert_dist(face[v1]));
2787  for (int d = 0; d < 3; d++)
2788  {
2789  pts[j][d] = (1-t)*verts(d,face[v0]) + t*verts(d,face[v1]);
2790  }
2791  if (sc != 0.0)
2792  {
2793  pts[j][3] = (1-t)*vals(face[v0]) + t*vals(face[v1]);
2794  }
2795  }
2796  if (sc != 0.0)
2797  {
2798  const double *dir = CuttingPlane->Equation();
2799  for (int j = 0; j < n; j++)
2800  {
2801  // dir points into the visible side, so we add a minus to val:
2802  const double val = -sc * (pts[j][3] - minv) / (maxv - minv);
2803  for (int d = 0; d < 3; d++)
2804  {
2805  pts[j][d] += val*dir[d];
2806  }
2807  }
2808  }
2809 
2810  target.addLine(gl3::Vertex::create(pts[0]), gl3::Vertex::create(pts[1]));
2811  if (n == 4)
2812  {
2813  target.addLine(gl3::Vertex::create(pts[2]), gl3::Vertex::create(pts[3]));
2814  }
2815  }
2816 }
2817 
2818 void VisualizationSceneSolution3d::PrepareCuttingPlane()
2819 {
2820  cplane_buf.clear();
2821  if (cp_drawelems && cplane && mesh->Dimension() == 3)
2822  {
2823  if (cplane == 2)
2824  {
2825  PrepareCuttingPlane2();
2826  }
2827  else if (cp_algo == 1)
2828  {
2829  CuttingPlaneFunc(1);
2830  }
2831  else
2832  {
2833  Vector vals, vert_dist;
2834  DenseMatrix pointmat;
2835  for (int i = 0; i < mesh->GetNE(); i++)
2836  {
2837  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
2838  RefinedGeometry *RefG =
2839  GLVisGeometryRefiner.Refine(geom, TimesToRefine);
2840  GridF->GetValues(i, RefG->RefPts, vals, pointmat);
2841  vert_dist.SetSize(pointmat.Width());
2842  for (int j = 0; j < pointmat.Width(); j++)
2843  {
2844  vert_dist(j) = CuttingPlane->Transform(&pointmat(0,j));
2845  }
2846  Array<int> &RG = RefG->RefGeoms;
2847 
2848  const int func = 0; // draw surface
2849  CutRefinedElement(cplane_buf, pointmat, vert_dist, vals, geom,
2850  RG, RG.Size()/Geometry::NumVerts[geom], func);
2851  }
2852  }
2853  }
2854  updated_bufs.emplace_back(&cplane_buf);
2855 }
2856 
2857 void VisualizationSceneSolution3d::PrepareCuttingPlane2()
2858 {
2859  int i, j, n = 0;
2860  double p[4][3], c[4], *coord;
2861  DenseMatrix pointmat, normals;
2862  Vector values;
2863  RefinedGeometry *RefG;
2864 
2865  Array<int> nodes;
2866  Array<int> partition (mesh -> GetNE());
2867 
2868  for (i = 0; i < mesh -> GetNE(); i++)
2869  {
2870  n = 0; // n will be the number of nodes behind the cutting plane
2871  mesh -> GetElementVertices(i, nodes);
2872  for (j = 0; j < nodes.Size(); j++)
2873  {
2874  if (node_pos[nodes[j]] >= 0.0) { n++; }
2875  }
2876  partition[i] = (n == nodes.Size()) ? 0 : 1;
2877  }
2878 
2879  for (i = 0; i < mesh -> GetNFaces(); i++)
2880  {
2881  int e1, e2;
2882  mesh -> GetFaceElements (i, &e1, &e2);
2883  if (e2 >= 0 && partition[e1] != partition[e2])
2884  {
2885  if (shading != 2)
2886  {
2887  mesh -> GetFaceVertices (i, nodes);
2888  for (j = 0; j < nodes.Size(); j++)
2889  {
2890  coord = mesh -> GetVertex(nodes[j]);
2891  p[j][0] = coord[0];
2892  p[j][1] = coord[1];
2893  p[j][2] = coord[2];
2894  c[j] = (*sol)(nodes[j]);
2895  }
2896 
2897  if (nodes.Size() == 3)
2898  {
2899  DrawTriangle(cplane_buf, p, c, minv, maxv);
2900  }
2901  else
2902  {
2903  DrawQuad(cplane_buf, p, c, minv, maxv);
2904  }
2905  }
2906  else // shading == 2
2907  {
2908  RefG = GLVisGeometryRefiner.Refine(mesh -> GetFaceBaseGeometry (i),
2909  TimesToRefine);
2910  // partition[e1] is 0 if e1 is behind the cutting plane
2911  // and 1 otherwise
2912  int dir = partition[e1];
2913  GridF -> GetFaceValues (i, dir, RefG->RefPts, values, pointmat);
2914  GetFaceNormals(i, dir, RefG->RefPts, normals);
2915  switch (mesh -> GetFaceBaseGeometry (i))
2916  {
2917  case Geometry::TRIANGLE: n = 3; break;
2918  case Geometry::SQUARE: n = 4; break;
2919  default:
2920  MFEM_ABORT("Invalid element type");
2921  break;
2922  }
2923  // DrawRefinedSurf (n, pointmat, values, RefG->RefGeoms);
2924  DrawPatch(cplane_buf, pointmat, values, normals, n, RefG->RefGeoms,
2925  minv, maxv, dir ? -3 : 2);
2926  } // end shading == 2
2927  }
2928  }
2929 }
2930 
2931 void VisualizationSceneSolution3d::PrepareCuttingPlaneLines()
2932 {
2933  cplines_buf.clear();
2934 
2935  if (cp_drawmesh && cplane && mesh->Dimension() == 3)
2936  {
2937  if (cplane == 2 && cp_drawmesh != 3)
2938  {
2939  PrepareCuttingPlaneLines2();
2940  }
2941  else
2942  {
2943  if (cp_drawmesh == 1 && cp_algo == 1)
2944  {
2945  CuttingPlaneFunc(2);
2946  }
2947  else if (cp_drawmesh == 1)
2948  {
2949  Vector vert_dist, vals;
2950  DenseMatrix pointmat;
2951  int num_faces = mesh->GetNFaces();
2952  if (mesh->NURBSext)
2953  {
2954  // Note: for NURBS meshes, the methods
2955  // Mesh::GetFaceTransformation() and
2956  // GridFunction::GetFaceValues() are not supported.
2957  cout << _MFEM_FUNC_NAME
2958  << ": NURBS mesh: cut faces will not be drawn!" << endl;
2959  num_faces = 0;
2960  }
2961  for (int i = 0; i < num_faces; i++)
2962  {
2963  const Geometry::Type geom = mesh->GetFaceBaseGeometry(i);
2964  RefinedGeometry *RefG =
2965  GLVisGeometryRefiner.Refine(geom, TimesToRefine);
2966  if (FaceShiftScale == 0.0)
2967  {
2968  ElementTransformation *T = mesh->GetFaceTransformation(i);
2969  T->Transform(RefG->RefPts, pointmat);
2970  }
2971  else
2972  {
2973  const int side = 2;
2974  GridF->GetFaceValues(i, side, RefG->RefPts, vals, pointmat);
2975  // For discontinuous grid function, we should draw two edges.
2976  }
2977  vert_dist.SetSize(pointmat.Width());
2978  for (int j = 0; j < pointmat.Width(); j++)
2979  {
2980  vert_dist(j) = CuttingPlane->Transform(&pointmat(0,j));
2981  }
2982  Array<int> &RG = RefG->RefGeoms;
2983 
2984  CutRefinedFace(cplines_buf, pointmat, vert_dist, vals, geom,
2985  RG, RG.Size()/Geometry::NumVerts[geom]);
2986  }
2987  }
2988  else if (cp_algo == 1)
2989  {
2990  CuttingPlaneFunc(3);
2991  }
2992  else
2993  {
2994  Vector vals, vert_dist;
2995  DenseMatrix pointmat;
2996  for (int i = 0; i < mesh->GetNE(); i++)
2997  {
2998  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
2999  RefinedGeometry *RefG =
3000  GLVisGeometryRefiner.Refine(geom, TimesToRefine);
3001  GridF->GetValues(i, RefG->RefPts, vals, pointmat);
3002  vert_dist.SetSize(pointmat.Width());
3003  for (int j = 0; j < pointmat.Width(); j++)
3004  {
3005  vert_dist(j) = CuttingPlane->Transform(&pointmat(0,j));
3006  }
3007  Array<int> &RG = RefG->RefGeoms;
3008 
3009  const int func = 1; // draw level lines
3010  CutRefinedElement(cplines_buf, pointmat, vert_dist, vals, geom,
3011  RG, RG.Size()/Geometry::NumVerts[geom], func);
3012  }
3013  }
3014  }
3015  }
3016 
3017  updated_bufs.emplace_back(&cplines_buf);
3018 }
3019 
3020 void VisualizationSceneSolution3d::PrepareCuttingPlaneLines2()
3021 {
3022  int i, j, n = 0;
3023  double point[4][4], *coord;
3024  DenseMatrix pointmat;
3025  Vector values;
3026  RefinedGeometry *RefG;
3027 
3028  Array<int> nodes;
3029  Array<int> partition (mesh -> GetNE());
3030  for (i = 0; i < mesh -> GetNE(); i++)
3031  {
3032  n = 0; // n will be the number of nodes behind the cutting plane
3033  mesh -> GetElementVertices(i,nodes);
3034  for (j=0; j<nodes.Size(); j++)
3035  {
3036  if (node_pos[nodes[j]] >= 0.0) { n++; }
3037  }
3038  partition[i] = (n == nodes.Size()) ? 0 : 1;
3039  }
3040 
3041  for (i = 0; i < mesh -> GetNFaces(); i++)
3042  {
3043  int e1, e2;
3044  mesh -> GetFaceElements (i, &e1, &e2);
3045  if (e2 >= 0 && partition[e1] != partition[e2])
3046  {
3047  if (shading != 2)
3048  {
3049  mesh -> GetFaceVertices (i, nodes);
3050  for (j = 0; j < nodes.Size(); j++)
3051  {
3052  coord = mesh -> GetVertex(nodes[j]);
3053  point[j][0] = coord[0];
3054  point[j][1] = coord[1];
3055  point[j][2] = coord[2];
3056  point[j][3] = (*sol)(nodes[j]);
3057  }
3058  gl3::GlBuilder line = cplines_buf.createBuilder();
3059  switch (cp_drawmesh)
3060  {
3061  case 1:
3062  // glBegin(GL_POLYGON);
3063  line.glBegin(GL_LINE_LOOP);
3064  for (j = 0; j < nodes.Size(); j++)
3065  {
3066  line.glVertex3dv (point[j]);
3067  }
3068  line.glEnd();
3069  break;
3070  case 2:
3071  DrawPolygonLevelLines(line, point[0], nodes.Size(), level, false);
3072  break;
3073  }
3074  }
3075  else // shading == 2
3076  {
3077  RefG = GLVisGeometryRefiner.Refine(mesh -> GetFaceBaseGeometry (i),
3078  TimesToRefine);
3079  // partition[e1] is 0 if e1 is behind the cutting plane
3080  // and 1 otherwise
3081  int di = partition[e1];
3082  GridF -> GetFaceValues (i, di, RefG->RefPts, values, pointmat);
3083  switch (mesh -> GetFaceBaseGeometry (i))
3084  {
3085  case Geometry::TRIANGLE: n = 3; break;
3086  case Geometry::SQUARE: n = 4; break;
3087  default:
3088  MFEM_ABORT("Invalid element type");
3089  break;
3090  }
3091  switch (cp_drawmesh)
3092  {
3093  case 1:
3094  DrawRefinedSurfEdges (n, pointmat, values, RefG->RefEdges);
3095  break;
3096  case 2:
3097  DrawRefinedSurfLevelLines (n, pointmat, values,
3098  RefG->RefGeoms);
3099  break;
3100  }
3101  } // end shading == 2
3102  }
3103  }
3104 }
3105 
3106 thread_local int triangle_counter;
3107 thread_local int quad_counter;
3108 
3109 void VisualizationSceneSolution3d::DrawTetLevelSurf(
3110  gl3::GlDrawable& target,
3111  const DenseMatrix &verts, const Vector &vals, const int *ind,
3112  const Array<double> &surf_levels, const DenseMatrix *grad)
3113 {
3114  double t, lvl, normal[3], vert[4][3], norm[4][3];
3115  int i, j, l, pos[4];
3116  bool flipped;
3117 
3118  gl3::GlBuilder draw = target.createBuilder();
3119 
3120  for (l = 0; l < surf_levels.Size(); l++)
3121  {
3122  lvl = surf_levels[l];
3123 
3124  for (i = 0; i < 4; i++)
3125  {
3126  pos[i] = ind[i];
3127  }
3128  i = 0;
3129  j = 4;
3130  flipped = false;
3131  do
3132  {
3133  // (i < j) is true
3134  while (vals(pos[i]) < lvl)
3135  {
3136  i++;
3137  if (i == j)
3138  {
3139  goto step_one;
3140  }
3141  }
3142  // (i < j) && (vals[pos[i]] >= lvl) is true
3143  do
3144  {
3145  j--;
3146  if (i == j)
3147  {
3148  goto step_one;
3149  }
3150  }
3151  while (vals(pos[j]) >= lvl);
3152  // (i < j) && (vals[pos[i]] >= lvl) && (vals[pos[j]] < lvl) is true
3153  Swap<int>(pos[i], pos[j]);
3154  flipped = !flipped;
3155  i++;
3156  }
3157  while (i < j);
3158  step_one:
3159  if (flipped)
3160  {
3161  if (i >= 2)
3162  {
3163  Swap<int>(pos[0], pos[1]);
3164  }
3165  else
3166  {
3167  Swap<int>(pos[2], pos[3]);
3168  }
3169  }
3170 
3171  if (j == 3)
3172  {
3173  Swap<int>(pos[0], pos[3]);
3174  j = 1;
3175  }
3176  if (j == 1)
3177  {
3178  for (int k = 0; k < 3; k++)
3179  {
3180  int p0 = pos[0];
3181  int p1 = pos[k+1];
3182  t = (lvl - vals(p0)) / (vals(p1) - vals(p0));
3183  for (int d = 0; d < 3; d++)
3184  {
3185  vert[k][d] = (1.0 - t) * verts(d, p0) + t * verts(d, p1);
3186  }
3187  if (grad)
3188  for (int d = 0; d < 3; d++)
3189  {
3190  norm[k][d] = (1.0 - t) * (*grad)(d, p0) + t * (*grad)(d, p1);
3191  }
3192  }
3193 
3194  if (grad == NULL)
3195  {
3196  if (!Compute3DUnitNormal(vert[0], vert[1], vert[2], normal))
3197  {
3198  MySetColor(draw, lvl, minv, maxv);
3199  draw.glNormal3dv(normal);
3200  draw.glBegin(GL_TRIANGLES);
3201  for (int k = 0; k < 3; k++)
3202  {
3203  draw.glVertex3dv(vert[k]);
3204  }
3205  draw.glEnd();
3206  triangle_counter++;
3207  }
3208  }
3209  else
3210  {
3211  MySetColor(draw, lvl, minv, maxv);
3212  draw.glBegin(GL_TRIANGLES);
3213  for (int k = 0; k < 3; k++)
3214  {
3215  Normalize(norm[k]);
3216  draw.glNormal3dv(norm[k]);
3217  draw.glVertex3dv(vert[k]);
3218  }
3219  draw.glEnd();
3220  triangle_counter++;
3221  }
3222  }
3223  else if (j == 2)
3224  {
3225  static const int idx[4][2] = { {0, 2}, {0, 3}, {1, 3}, {1, 2} };
3226 
3227  for (int k = 0; k < 4; k++)
3228  {
3229  int p0 = pos[idx[k][0]];
3230  int p1 = pos[idx[k][1]];
3231  t = (lvl - vals(p0)) / (vals(p1) - vals(p0));
3232  for (int d = 0; d < 3; d++)
3233  {
3234  vert[k][d] = (1.0 - t) * verts(d, p0) + t * verts(d, p1);
3235  }
3236  if (grad)
3237  for (int d = 0; d < 3; d++)
3238  {
3239  norm[k][d] = (1.0 - t) * (*grad)(d, p0) + t * (*grad)(d, p1);
3240  }
3241  }
3242 
3243  if (grad == NULL)
3244  {
3245  if (!Compute3DUnitNormal(vert[0], vert[1], vert[2], vert[3],
3246  normal))
3247  {
3248  MySetColor(draw, lvl, minv, maxv);
3249  draw.glNormal3dv(normal);
3250  draw.glBegin(GL_QUADS);
3251  for (int k = 0; k < 4; k++)
3252  {
3253  draw.glVertex3dv(vert[k]);
3254  }
3255  draw.glEnd();
3256  quad_counter++;
3257  }
3258  }
3259  else
3260  {
3261  MySetColor(draw, lvl, minv, maxv);
3262  draw.glBegin(GL_QUADS);
3263  for (int k = 0; k < 4; k++)
3264  {
3265  Normalize(norm[k]);
3266  draw.glNormal3dv(norm[k]);
3267  draw.glVertex3dv(vert[k]);
3268  }
3269  draw.glEnd();
3270  quad_counter++;
3271  }
3272  }
3273  }
3274 }
3275 
3276 // static method
3277 int VisualizationSceneSolution3d::GetPyramidFaceSplits(
3278  const Array<bool> &quad_diag, const Array<int> &faces,
3279  const Array<int> &ofaces)
3280 {
3281  int fs = 0;
3282  bool diag = quad_diag[faces[0]];
3283  if ((ofaces[0]/2)%2) // orientations 2,3,6,7
3284  {
3285  diag = !diag;
3286  }
3287  fs = 2*fs + diag;
3288  return fs;
3289 }
3290 
3291 void VisualizationSceneSolution3d::DrawRefinedPyramidLevelSurf(
3292  gl3::GlDrawable& target,
3293  const DenseMatrix &verts, const Vector &vals, const int *RG, const int np,
3294  const int face_splits, const DenseMatrix *grad)
3295 {
3296 #if 1
3297  static const int pyr_tets[2][4] =
3298  {
3299  { 0, 1, 2, 4 }, { 0, 2, 3, 4 }
3300  };
3301  for (int k = 0; k < np; k++)
3302  {
3303  const int *hv = &RG[5*k];
3304  if (hv[4] > 0)
3305  {
3306  for (int j = 0; j < 2; j++)
3307  {
3308  int m_ind[4];
3309  for (int i = 0; i < 4; i++)
3310  {
3311  m_ind[i] = hv[pyr_tets[j][i]];
3312  }
3313  DrawTetLevelSurf(target, verts, vals, m_ind, levels, grad);
3314  }
3315  }
3316  else
3317  {
3318  DrawTetLevelSurf(target, verts, vals, hv, levels, grad);
3319  }
3320  }
3321 #else
3322  // MLS: Not sure how to adapt this to Pyramids so skip it for now
3323  static const int pri_tets[8-2][3][4] =
3324  {
3325  // 0 = 000 (see below; prism is split into 6 tets)
3326  { { 0, 1, 2, 5 }, { 0, 1, 5, 4 }, { 0, 3, 4, 5 } }, // 1 = 001
3327  { { 0, 1, 2, 4 }, { 0, 2, 3, 4 }, { 2, 3, 4, 5 } }, // 2 = 010
3328  { { 0, 1, 2, 4 }, { 0, 2, 5, 4 }, { 0, 3, 4, 5 } }, // 3 = 011
3329  { { 0, 1, 2, 3 }, { 1, 2, 3, 5 }, { 1, 3, 4, 5 } }, // 4 = 100
3330  { { 0, 1, 2, 5 }, { 0, 1, 5, 3 }, { 1, 3, 4, 5 } }, // 5 = 101
3331  { { 0, 1, 2, 3 }, { 1, 2, 3, 4 }, { 2, 3, 4, 5 } } // 6 = 110
3332  // 7 = 111 (see below; prism is split into 6 tets)
3333  };
3334  static const int pri_tets_0[6][4] =
3335  { {0,1,2,6}, {0,1,6,4}, {0,2,3,6}, {0,6,3,4}, {2,3,6,5}, {3,4,6,5} };
3336  static const int pri_tets_7[6][4] =
3337  { {0,1,2,6}, {0,2,5,6}, {0,1,6,3}, {0,3,6,5}, {1,3,4,6}, {3,4,6,5} };
3338  const int fs2 = (~face_splits&7)/2 + 4*((~face_splits&7)%2);
3339  const int n = (np == 1) ? 1 : TimesToRefine;
3340 
3341  double vs_data[7], pm_data[3*7], gd_data[3*7];
3342  Vector vs(vs_data, 7);
3343  DenseMatrix pm(pm_data, 3, 7), gd(gd_data, 3, 7);
3344 
3345  for (int k = 0, l0 = -1, l1 = -1; k < np; k++)
3346  {
3347  const int pk = k % (n*n);
3348  if (pk == 0)
3349  {
3350  l0 = 0; l1 = 2*n-1;
3351  }
3352  else if (pk == l1)
3353  {
3354  const int s = l1-l0;
3355  l0 = l1;
3356  l1 += (s-2);
3357  }
3358  const int fsl = ((pk-l0)%2 == 0) ? face_splits : fs2;
3359  // The algorithm for choosing 'fsl' used above assumes the refined prisms
3360  // are listed in certain order -- see the prism case in
3361  // mfem::GeometryRefiner::Refine.
3362  const int *pv = &RG[6*k];
3363  if (fsl == 0)
3364  {
3365  for (int j = 0; j < 6; j++)
3366  {
3367  const int idx = pv[j];
3368  for (int d = 0; d < 3; d++)
3369  {
3370  pm(d,j) = verts(d,idx);
3371  if (grad) { gd(d,j) = (*grad)(d,idx); }
3372  }
3373  vs(j) = vals(idx);
3374  }
3375  for (int d = 0; d < 3; d++)
3376  {
3377  pm(d,6) = 0.5*(pm(d,1) + pm(d,5));
3378  if (grad) { gd(d,6) = 0.5*(gd(d,1) + gd(d,5)); }
3379  }
3380  vs(6) = 0.5*(vs(1) + vs(5));
3381  const DenseMatrix *gd_ = grad ? &gd : NULL;
3382  for (int k = 0; k < 6; k++)
3383  {
3384  DrawTetLevelSurf(pm, vs, pri_tets_0[k], levels, gd_);
3385  }
3386  }
3387  else if (fsl == 7)
3388  {
3389  for (int j = 0; j < 6; j++)
3390  {
3391  const int idx = pv[j];
3392  for (int d = 0; d < 3; d++)
3393  {
3394  pm(d,j) = verts(d,idx);
3395  if (grad) { gd(d,j) = (*grad)(d,idx); }
3396  }
3397  vs(j) = vals(idx);
3398  }
3399  for (int d = 0; d < 3; d++)
3400  {
3401  pm(d,6) = 0.5*(pm(d,2) + pm(d,4));
3402  if (grad) { gd(d,6) = 0.5*(gd(d,2) + gd(d,4)); }
3403  }
3404  vs(6) = 0.5*(vs(2) + vs(4));
3405  const DenseMatrix *gd_ = grad ? &gd : NULL;
3406  for (int k = 0; k < 6; k++)
3407  {
3408  DrawTetLevelSurf(pm, vs, pri_tets_7[k], levels, gd_);
3409  }
3410  }
3411  else
3412  {
3413  int m_ind[4];
3414  for (int j = 0; j < 3; j++)
3415  {
3416  for (int i = 0; i < 4; i++)
3417  {
3418  m_ind[i] = pv[pri_tets[fsl-1][j][i]];
3419  }
3420  DrawTetLevelSurf(verts, vals, m_ind, levels, grad);
3421  }
3422  }
3423  }
3424 #endif
3425 }
3426 
3427 // static method
3428 int VisualizationSceneSolution3d::GetWedgeFaceSplits(
3429  const Array<bool> &quad_diag, const Array<int> &faces,
3430  const Array<int> &ofaces)
3431 {
3432  int fs = 0;
3433  for (int lf = 2; lf < 5; lf++)
3434  {
3435  bool diag = quad_diag[faces[lf]];
3436  if ((ofaces[lf]/2)%2) // orientations 2,3,6,7
3437  {
3438  diag = !diag;
3439  }
3440  fs = 2*fs + diag;
3441  }
3442  return fs;
3443 }
3444 
3445 void VisualizationSceneSolution3d::DrawRefinedWedgeLevelSurf(
3446  gl3::GlDrawable& target,
3447  const DenseMatrix &verts, const Vector &vals, const int *RG, const int np,
3448  const int face_splits, const DenseMatrix *grad)
3449 {
3450 #if 0
3451  static const int pri_tets[3][4] =
3452  {
3453  { 0, 1, 2, 5 }, { 0, 1, 5, 3 }, { 1, 3, 4, 5 }
3454  };
3455  for (int k = 0; k < np; k++)
3456  {
3457  const int *hv = &RG[6*k];
3458  for (int j = 0; j < 3; j++)
3459  {
3460  int m_ind[4];
3461  for (int i = 0; i < 4; i++)
3462  {
3463  m_ind[i] = hv[pri_tets[j][i]];
3464  }
3465  DrawTetLevelSurf(verts, vals, m_ind, levels, grad);
3466  }
3467  }
3468 #else
3469  static const int pri_tets[8-2][3][4] =
3470  {
3471  // 0 = 000 (see below; prism is split into 6 tets)
3472  { { 0, 1, 2, 5 }, { 0, 1, 5, 4 }, { 0, 3, 4, 5 } }, // 1 = 001
3473  { { 0, 1, 2, 4 }, { 0, 2, 3, 4 }, { 2, 3, 4, 5 } }, // 2 = 010
3474  { { 0, 1, 2, 4 }, { 0, 2, 5, 4 }, { 0, 3, 4, 5 } }, // 3 = 011
3475  { { 0, 1, 2, 3 }, { 1, 2, 3, 5 }, { 1, 3, 4, 5 } }, // 4 = 100
3476  { { 0, 1, 2, 5 }, { 0, 1, 5, 3 }, { 1, 3, 4, 5 } }, // 5 = 101
3477  { { 0, 1, 2, 3 }, { 1, 2, 3, 4 }, { 2, 3, 4, 5 } } // 6 = 110
3478  // 7 = 111 (see below; prism is split into 6 tets)
3479  };
3480  static const int pri_tets_0[6][4] =
3481  { {0,1,2,6}, {0,1,6,4}, {0,2,3,6}, {0,6,3,4}, {2,3,6,5}, {3,4,6,5} };
3482  static const int pri_tets_7[6][4] =
3483  { {0,1,2,6}, {0,2,5,6}, {0,1,6,3}, {0,3,6,5}, {1,3,4,6}, {3,4,6,5} };
3484  const int fs2 = (~face_splits&7)/2 + 4*((~face_splits&7)%2);
3485  const int n = (np == 1) ? 1 : TimesToRefine;
3486 
3487  double vs_data[7], pm_data[3*7], gd_data[3*7];
3488  Vector vs(vs_data, 7);
3489  DenseMatrix pm(pm_data, 3, 7), gd(gd_data, 3, 7);
3490 
3491  for (int k = 0, l0 = -1, l1 = -1; k < np; k++)
3492  {
3493  const int pk = k % (n*n);
3494  if (pk == 0)
3495  {
3496  l0 = 0; l1 = 2*n-1;
3497  }
3498  else if (pk == l1)
3499  {
3500  const int s = l1-l0;
3501  l0 = l1;
3502  l1 += (s-2);
3503  }
3504  const int fsl = ((pk-l0)%2 == 0) ? face_splits : fs2;
3505  // The algorithm for choosing 'fsl' used above assumes the refined prisms
3506  // are listed in certain order -- see the prism case in
3507  // mfem::GeometryRefiner::Refine.
3508  const int *pv = &RG[6*k];
3509  if (fsl == 0)
3510  {
3511  for (int j = 0; j < 6; j++)
3512  {
3513  const int idx = pv[j];
3514  for (int d = 0; d < 3; d++)
3515  {
3516  pm(d,j) = verts(d,idx);
3517  if (grad) { gd(d,j) = (*grad)(d,idx); }
3518  }
3519  vs(j) = vals(idx);
3520  }
3521  for (int d = 0; d < 3; d++)
3522  {
3523  pm(d,6) = 0.5*(pm(d,1) + pm(d,5));
3524  if (grad) { gd(d,6) = 0.5*(gd(d,1) + gd(d,5)); }
3525  }
3526  vs(6) = 0.5*(vs(1) + vs(5));
3527  const DenseMatrix *gd_ = grad ? &gd : NULL;
3528  for (int j = 0; j < 6; j++)
3529  {
3530  DrawTetLevelSurf(target, pm, vs, pri_tets_0[j], levels, gd_);
3531  }
3532  }
3533  else if (fsl == 7)
3534  {
3535  for (int j = 0; j < 6; j++)
3536  {
3537  const int idx = pv[j];
3538  for (int d = 0; d < 3; d++)
3539  {
3540  pm(d,j) = verts(d,idx);
3541  if (grad) { gd(d,j) = (*grad)(d,idx); }
3542  }
3543  vs(j) = vals(idx);
3544  }
3545  for (int d = 0; d < 3; d++)
3546  {
3547  pm(d,6) = 0.5*(pm(d,2) + pm(d,4));
3548  if (grad) { gd(d,6) = 0.5*(gd(d,2) + gd(d,4)); }
3549  }
3550  vs(6) = 0.5*(vs(2) + vs(4));
3551  const DenseMatrix *gd_ = grad ? &gd : NULL;
3552  for (int j = 0; j < 6; j++)
3553  {
3554  DrawTetLevelSurf(target, pm, vs, pri_tets_7[j], levels, gd_);
3555  }
3556  }
3557  else
3558  {
3559  int m_ind[4];
3560  for (int j = 0; j < 3; j++)
3561  {
3562  for (int i = 0; i < 4; i++)
3563  {
3564  m_ind[i] = pv[pri_tets[fsl-1][j][i]];
3565  }
3566  DrawTetLevelSurf(target, verts, vals, m_ind, levels, grad);
3567  }
3568  }
3569  }
3570 #endif
3571 }
3572 
3573 // static method
3574 int VisualizationSceneSolution3d::GetHexFaceSplits(
3575  const Array<bool> &quad_diag, const Array<int> &faces,
3576  const Array<int> &ofaces)
3577 {
3578  int fs = 0;
3579  for (int lf = 0; lf < 6; lf++)
3580  {
3581  bool diag = quad_diag[faces[lf]];
3582  if ((ofaces[lf]/2)%2) // orientations 2,3,6,7
3583  {
3584  diag = !diag;
3585  }
3586  fs = 2*fs + diag;
3587  }
3588  return fs;
3589 }
3590 
3591 void VisualizationSceneSolution3d::DrawRefinedHexLevelSurf(
3592  gl3::GlDrawable& target,
3593  const DenseMatrix &verts, const Vector &vals, const int *RG, const int nh,
3594  const int face_splits, const DenseMatrix *grad)
3595 {
3596 #if 0
3597  static const int hex_tets[6][4] =
3598  {
3599  { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
3600  { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
3601  };
3602  for (int k = 0; k < nh; k++)
3603  {
3604  const int *hv = &RG[8*k];
3605  for (int j = 0; j < 6; j++)
3606  {
3607  int m_ind[4];
3608  for (int i = 0; i < 4; i++)
3609  {
3610  m_ind[i] = hv[hex_tets[j][i]];
3611  }
3612  DrawTetLevelSurf(verts, vals, m_ind, levels, grad);
3613  }
3614  }
3615 #else
3616  const int n = (nh == 1) ? 1 : TimesToRefine;
3617 
3618  double vs_data[9], pm_data[3*9], gd_data[3*9];
3619  Vector vs(vs_data, 9);
3620  DenseMatrix pm(pm_data, 3, 9), gd(gd_data, 3, 9);
3621 
3622  for (int k = 0; k < nh; k++)
3623  {
3624  const int ix = k%n, iy = (k/n)%n, iz = k/(n*n);
3625  int fsl = face_splits;
3626  if (ix != 0)
3627  {
3628  // Copy the right face bit (face 2, bit mask 8) to the left face bit
3629  // (face 4, bit mask 2) and flip it.
3630  fsl = (fsl&(63-2)) + (2-(fsl&8)/4);
3631  }
3632  if (iy != 0)
3633  {
3634  // Copy the back face bit (face 3, bit mask 4) to the front face bit
3635  // (face 1, bit mask 16) and flip it.
3636  fsl = (fsl&(63-16)) + (16-(fsl&4)*4);
3637  }
3638  if (iz != 0)
3639  {
3640  // Copy the top face bit (face 5, bit mask 1) to the bottom face bit
3641  // (face 0, bit mask 32) and flip it.
3642  fsl = (fsl&(63-32)) + (32-(fsl&1)*32);
3643  }
3644 
3645  const int *hv_orig = &RG[8*k], *hv;
3646  int hv1[8], hv2[8];
3647 #if 1
3648  // Find a pair of opposite faces that are split into triangles using
3649  // two parallel diagonals (if such pair exists).
3650  if (1-(fsl&1) == (fsl&32)/32)
3651  {
3652  // The bottom and top faces are split in the same direction.
3653  hv = hv_orig;
3654  }
3655  else if (4-(fsl&4) == (fsl&16)/4)
3656  {
3657  // The front and back faces are split in the same direction.
3658  // Rotate the hex around the x-axis to bring front-back to bottom-top.
3659  static const int rot_x[8] = { 4, 5, 1, 0, 7, 6, 2, 3 };
3660  for (int j = 0; j < 8; j++)
3661  {
3662  hv1[j] = hv_orig[rot_x[j]];
3663  }
3664  hv = hv1;
3665  // fsl bits change: a|b|c|d|e|f -> b|f|1-c|a|1-e|d
3666  fsl = (fsl&16)*2 + (fsl&1)*16 + (8-(fsl&8)) + (fsl&32)/8 +
3667  (2-(fsl&2)) + (fsl&4)/4;
3668  }
3669  else if (2-(fsl&2) == (fsl&8)/4)
3670  {
3671  // The left and right faces are split in the same direction.
3672  // Rotate the hex around the y-axis to bring left-right to top-bottom.
3673  static const int rot_y[8] = { 1, 5, 6, 2, 0, 4, 7, 3 };
3674  for (int j = 0; j < 8; j++)
3675  {
3676  hv1[j] = hv_orig[rot_y[j]];
3677  }
3678  hv = hv1;
3679  // fsl bit change: a|b|c|d|e|f -> 1-c|1-b|1-f|1-d|1-a|1-e
3680  fsl = ~fsl&63;
3681  fsl = (fsl&8)*4 + (fsl&16) + (fsl&1)*8 + (fsl&4) + (fsl&32)/16 +
3682  (fsl&2)/2;
3683  }
3684  else
3685 #endif
3686  {
3687  // All opposite faces are split in opposite directions.
3688  // Split the hex into 12 tets using the hex center as a vertex in all
3689  // 12 tets.
3690  for (int j = 0; j < 8; j++)
3691  {
3692  const int idx = hv_orig[j];
3693  for (int d = 0; d < 3; d++)
3694  {
3695  pm(d,j) = verts(d,idx);
3696  if (grad) { gd(d,j) = (*grad)(d,idx); }
3697  }
3698  vs(j) = vals(idx);
3699  }
3700  for (int d = 0; d < 3; d++)
3701  {
3702  pm(d,8) = 0.0;
3703  if (grad) { gd(d,8) = 0.0; }
3704  }
3705  vs(8) = 0.0;
3706  for (int j = 0; j < 8; j++)
3707  {
3708  for (int d = 0; d < 3; d++)
3709  {
3710  pm(d,8) += pm(d,j);
3711  if (grad) { gd(d,8) += gd(d,j); }
3712  }
3713  vs(8) += vs(j);
3714  }
3715  for (int d = 0; d < 3; d++)
3716  {
3717  pm(d,8) *= 0.125;
3718  if (grad) { gd(d,8) *= 0.125; }
3719  }
3720  vs(8) *= 0.125;
3721 
3722  typedef Mesh::hex_t hex_t;
3723  for (int j = 0; j < hex_t::NumFaces; j++)
3724  {
3725  int tv[8];
3726  const int *fv = hex_t::FaceVert[j];
3727  const bool diag = fsl&(32>>j);
3728  if (diag == 0)
3729  {
3730  tv[0] = fv[2];
3731  tv[1] = fv[1];
3732  tv[2] = fv[0];
3733  tv[3] = 8;
3734  tv[4] = fv[0];
3735  tv[5] = fv[3];
3736  tv[6] = fv[2];
3737  tv[7] = 8;
3738  }
3739  else
3740  {
3741  tv[0] = fv[1];
3742  tv[1] = fv[0];
3743  tv[2] = fv[3];
3744  tv[3] = 8;
3745  tv[4] = fv[3];
3746  tv[5] = fv[2];
3747  tv[6] = fv[1];
3748  tv[7] = 8;
3749  }
3750  const DenseMatrix *gp = grad ? &gd : NULL;
3751  DrawTetLevelSurf(target, pm, vs, &tv[0], levels, gp);
3752  DrawTetLevelSurf(target, pm, vs, &tv[4], levels, gp);
3753  }
3754 
3755  continue;
3756  }
3757 
3758  // Rotate the hex so that the diagonal edge splitting the bottom face is
3759  // the edge 0-2.
3760  if ((fsl&32) == 0)
3761  {
3762  // Rotate the hex around the z-axis -- left-right becomes front-back.
3763  static const int rot_z[8] = { 3, 0, 1, 2, 7, 4, 5, 6 };
3764  for (int j = 0; j < 8; j++)
3765  {
3766  hv2[j] = hv[rot_z[j]];
3767  }
3768  hv = hv2;
3769  // fsl bit change: a|b|c|d|e|f -> 1-a|e|b|c|d|1-f
3770  fsl = (32-(fsl&32)) + (fsl&2)*8 + (fsl&(16+8+4))/2 + (1-(fsl&1));
3771  }
3772 
3773  // Split the hex into two prisms using the diagonal face 0-2-6-4.
3774  const int pv[2][6] =
3775  {
3776  { hv[0], hv[1], hv[2], hv[4], hv[5], hv[6] },
3777  { hv[2], hv[3], hv[0], hv[6], hv[7], hv[4] }
3778  };
3779  // Choose the shorter diagonal on the face 0-2-6-4.
3780  const double l06 = Distance(&verts(0,hv[0]), &verts(0,hv[6]), 3);
3781  const double l24 = Distance(&verts(0,hv[2]), &verts(0,hv[4]), 3);
3782  const bool diag = (l06 > 1.01*l24);
3783  const int fs1 = (fsl&(16+8))/4 + !diag; // a|b|c|d|e|f -> b|c|1-diag
3784  const int fs2 = (fsl&(4+2)) + diag; // a|b|c|d|e|f -> d|e|diag
3785  DrawRefinedWedgeLevelSurf(target, verts, vals, pv[0], 1, fs1, grad);
3786  DrawRefinedWedgeLevelSurf(target, verts, vals, pv[1], 1, fs2, grad);
3787  }
3788 #endif
3789 }
3790 
3791 void VisualizationSceneSolution3d::PrepareLevelSurf()
3792 {
3793  static const int ident[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
3794 
3795  Vector vals;
3796  DenseMatrix pointmat, grad;
3797  Array<int> vertices;
3798 
3799  lsurf_buf.clear();
3800  if (drawlsurf == 0 || mesh->Dimension() != 3)
3801  {
3802  // Create empty list
3803  return;
3804  }
3805 
3806  triangle_counter = quad_counter = 0;
3807 
3808  levels.SetSize(nlevels);
3809  for (int l = 0; l < nlevels; l++)
3810  {
3811  double lvl = ((double)(50*l+drawlsurf) / (nlevels*50));
3812  levels[l] = ULogVal(lvl);
3813  }
3814 
3815  // For every quad face, choose the shorter diagonal to split the quad into
3816  // two triangles. Elements adjacent to that quad face (wedge or hex) will use
3817  // the same diagonal when subdividing the element.
3818  Array<bool> quad_diag;
3819  if (mesh->HasGeometry(Geometry::SQUARE))
3820  {
3821  quad_diag.SetSize(mesh->GetNFaces());
3822  for (int fi = 0; fi < mesh->GetNFaces(); fi++)
3823  {
3824  const Element *face = mesh->GetFace(fi);
3825  if (face->GetType() != Element::QUADRILATERAL) { continue; }
3826  ElementTransformation *T = mesh->GetFaceTransformation(fi);
3827  T->Transform(*Geometries.GetVertices(Geometry::SQUARE), pointmat);
3828  const double l02 = Distance(&pointmat(0,0), &pointmat(0,2), 3);
3829  const double l13 = Distance(&pointmat(0,1), &pointmat(0,3), 3);
3830  quad_diag[fi] = (l02 > 1.01*l13);
3831  }
3832  }
3833 
3834  Array<int> faces, ofaces;
3835 
3836  if (shading != 2)
3837  {
3838  for (int ie = 0; ie < mesh->GetNE(); ie++)
3839  {
3840  mesh->GetPointMatrix(ie, pointmat);
3841  mesh->GetElementVertices(ie, vertices);
3842  vals.SetSize(vertices.Size());
3843  for (int j = 0; j < vertices.Size(); j++)
3844  {
3845  vals(j) = (*sol)(vertices[j]);
3846  }
3847 
3848  switch (mesh->GetElementType(ie))
3849  {
3850  case Element::TETRAHEDRON:
3851  DrawTetLevelSurf(lsurf_buf, pointmat, vals, ident, levels);
3852  break;
3853  case Element::PYRAMID:
3854  {
3855  mesh->GetElementFaces(ie, faces, ofaces);
3856  const int fs = GetPyramidFaceSplits(quad_diag, faces, ofaces);
3857  DrawRefinedPyramidLevelSurf(lsurf_buf, pointmat, vals, ident, 1, fs);
3858  }
3859  case Element::WEDGE:
3860  {
3861  mesh->GetElementFaces(ie, faces, ofaces);
3862  const int fs = GetWedgeFaceSplits(quad_diag, faces, ofaces);
3863  DrawRefinedWedgeLevelSurf(lsurf_buf, pointmat, vals, ident, 1, fs);
3864  }
3865  break;
3866  case Element::HEXAHEDRON:
3867  {
3868  mesh->GetElementFaces(ie, faces, ofaces);
3869  const int fs = GetHexFaceSplits(quad_diag, faces, ofaces);
3870  DrawRefinedHexLevelSurf(lsurf_buf, pointmat, vals, ident, 1, fs);
3871  }
3872  break;
3873  default:
3874  MFEM_ABORT("Unrecognized 3D element type \""
3875  << mesh->GetElementType(ie) << "\"");
3876  }
3877  }
3878  }
3879  else // shading == 2
3880  {
3881  RefinedGeometry *RefG;
3882 #define GLVIS_SMOOTH_LEVELSURF_NORMALS
3883 #ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS
3884  const DenseMatrix *gp = &grad;
3885 #else
3886  const DenseMatrix *gp = NULL;
3887 #endif
3888 
3889  for (int ie = 0; ie < mesh->GetNE(); ie++)
3890  {
3891  const Geometry::Type geom = mesh->GetElementBaseGeometry(ie);
3892 
3893  RefG = GLVisGeometryRefiner.Refine(geom, TimesToRefine);
3894  GridF->GetValues(ie, RefG->RefPts, vals, pointmat);
3895 #ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS
3896  GridF->GetGradients(ie, RefG->RefPts, grad);
3897 #endif
3898 
3899  Array<int> &RG = RefG->RefGeoms;
3900  const int nv = mesh->GetElement(ie)->GetNVertices();
3901  const int nre = RG.Size()/nv;
3902 
3903  if (geom == Geometry::TETRAHEDRON)
3904  {
3905  for (int k = 0; k < nre; k++)
3906  {
3907  DrawTetLevelSurf(lsurf_buf, pointmat, vals, &RG[nv*k], levels, gp);
3908  }
3909  }
3910  else if (geom == Geometry::PYRAMID)
3911  {
3912  mesh->GetElementFaces(ie, faces, ofaces);
3913  const int fs = GetPyramidFaceSplits(quad_diag, faces, ofaces);
3914  DrawRefinedPyramidLevelSurf(lsurf_buf, pointmat, vals, RG, nre, fs, gp);
3915  }
3916  else if (geom == Geometry::PRISM)
3917  {
3918  mesh->GetElementFaces(ie, faces, ofaces);
3919  const int fs = GetWedgeFaceSplits(quad_diag, faces, ofaces);
3920  DrawRefinedWedgeLevelSurf(lsurf_buf, pointmat, vals, RG, nre, fs, gp);
3921  }
3922  else if (geom == Geometry::CUBE)
3923  {
3924  mesh->GetElementFaces(ie, faces, ofaces);
3925  const int fs = GetHexFaceSplits(quad_diag, faces, ofaces);
3926  DrawRefinedHexLevelSurf(lsurf_buf, pointmat, vals, RG, nre, fs, gp);
3927  }
3928  }
3929  }
3930 
3931  updated_bufs.emplace_back(&lsurf_buf);
3932 
3933 #ifdef GLVIS_DEBUG
3934  cout << "VisualizationSceneSolution3d::PrepareLevelSurf() : "
3935  << triangle_counter << " triangles + " << quad_counter
3936  << " quads used" << endl;
3937 #endif
3938 }
3939 
3940 gl3::SceneInfo VisualizationSceneSolution3d::GetSceneObjs()
3941 {
3942  if (colorbar)
3943  {
3944  Array<double>* cb_level = nullptr,
3945  * cb_levels = nullptr;
3946  if (drawlsurf) { cb_levels = &levels; }
3947  if (drawmesh == 2 || cp_drawmesh >= 2) { cb_level = &level; }
3948  PrepareColorBar(minv, maxv, cb_level, cb_levels);
3949  }
3950  gl3::SceneInfo scene = VisualizationSceneScalarData::GetSceneObjs();
3951  gl3::RenderParams params = GetMeshDrawParams();
3952  params.use_clip_plane = cplane;
3953  double* cp_eqn = CuttingPlane->Equation();
3954  params.clip_plane_eqn = {cp_eqn[0], cp_eqn[1], cp_eqn[2], cp_eqn[3]};
3955  params.contains_translucent = matAlpha < 1.0;
3956 
3957  if (drawlsurf)
3958  {
3959  scene.queue.emplace_back(params, &lsurf_buf);
3960  }
3961  if (drawelems)
3962  {
3963  scene.queue.emplace_back(params, &disp_buf);
3964  }
3965  if (cplane && cp_drawelems)
3966  {
3967  params.use_clip_plane = false;
3968  scene.queue.emplace_back(params, &cplane_buf);
3969  params.use_clip_plane = true;
3970  }
3971  params.contains_translucent = false;
3972  // draw orderings -- color modes
3973  if (draworder == 1)
3974  {
3975  scene.queue.emplace_back(params, &order_noarrow_buf);
3976  }
3977  else if (draworder == 2)
3978  {
3979  scene.queue.emplace_back(params, &order_buf);
3980  }
3981  params.mesh_material = VisualizationScene::BLK_MAT;
3982  // everything below will be drawn in "black"
3983  params.static_color = GetLineColor();
3984  params.num_pt_lights = 0;
3985  if (drawmesh)
3986  {
3987  scene.queue.emplace_back(params, &line_buf);
3988  }
3989  if (cp_drawmesh)
3990  {
3991  params.use_clip_plane = false;
3992  scene.queue.emplace_back(params, &cplines_buf);
3993  }
3994  // draw orderings -- "black" modes
3995  if (draworder == 3)
3996  {
3997  scene.queue.emplace_back(params, &order_noarrow_buf);
3998  }
3999  else if (draworder == 4)
4000  {
4001  scene.queue.emplace_back(params, &order_buf);
4002  }
4003  return scene;
4004 }
4005 
4006 void VisualizationSceneSolution3d::glTF_Export()
4007 {
4008  string name = "GLVis_scene_000";
4009 
4010  glTF_Builder bld(name);
4011 
4012  auto palette_mat = AddPaletteMaterial(bld);
4013  auto black_mat = AddBlackMaterial(bld);
4014  auto buf = bld.addBuffer("buffer");
4015  if (drawelems) { glTF_ExportElements(bld, buf, palette_mat, disp_buf); }
4016  if (drawmesh) { glTF_ExportMesh(bld, buf, black_mat, line_buf); }
4017  if (cplane && cp_drawelems)
4018  {
4019  auto cp_elems_node = AddModelNode(bld, "CP Elements");
4020  auto cp_elems_mesh = bld.addMesh("CP Elements Mesh");
4021  bld.addNodeMesh(cp_elems_node, cp_elems_mesh);
4022 
4023  int ntria = AddTriangles(
4024  bld,
4025  cp_elems_mesh,
4026  buf,
4027  palette_mat,
4028  cplane_buf);
4029  if (ntria == 0)
4030  {
4031  cout << "glTF export: no cp elements found to export!" << endl;
4032  }
4033  }
4034  if (cp_drawmesh)
4035  {
4036  auto cp_lines_node = AddModelNode(bld, "CP Lines");
4037  auto cp_lines_mesh = bld.addMesh("CP Lines Mesh");
4038  bld.addNodeMesh(cp_lines_node, cp_lines_mesh);
4039 
4040  int nlines = AddLines(
4041  bld,
4042  cp_lines_mesh,
4043  buf,
4044  black_mat,
4045  cplines_buf);
4046  if (nlines == 0)
4047  {
4048  cout << "glTF export: no cp mesh/level lines found to export!" << endl;
4049  }
4050  }
4051  if (drawlsurf)
4052  {
4053  auto lsurf_node = AddModelNode(bld, "Level Surface");
4054  auto lsurf_mesh = bld.addMesh("Level Surface Mesh");
4055  bld.addNodeMesh(lsurf_node, lsurf_mesh);
4056 
4057  int ntria = AddTriangles(
4058  bld,
4059  lsurf_mesh,
4060  buf,
4061  palette_mat,
4062  lsurf_buf);
4063  if (ntria == 0)
4064  {
4065  cout << "glTF export: no level surface elements found to export!"
4066  << endl;
4067  }
4068  }
4069  if (drawaxes) { glTF_ExportBox(bld, buf, black_mat); }
4070  bld.writeFile();
4071 
4072  cout << "Exported glTF -> " << name << ".gltf" << endl;
4073 }
VisualizationScene::cut_lambda
double cut_lambda
Amount of face cutting with keys Ctrl-F3/F4 (0: no cut, 1: cut to edges)
Definition: openglvis.hpp:187
SendExposeEvent
void SendExposeEvent()
Send expose event. In our case MyReshape is executed and Draw after it.
Definition: aux_vis.cpp:346
VisualizationScene::cut_updated
bool cut_updated
Have the reference geometries been updated for the cut?
Definition: openglvis.hpp:189
quad_counter
thread_local int quad_counter
Definition: vssolution3d.cpp:3107
KeyvPressed
void KeyvPressed()
Definition: vsvector.cpp:141
GLVisGeometryRefiner
thread_local GeometryRefiner GLVisGeometryRefiner
Definition: glvis.cpp:71
VisualizationSceneSolution3d::GetWedgeFaceSplits
static int GetWedgeFaceSplits(const Array< bool > &quad_diag, const Array< int > &faces, const Array< int > &ofaces)
Definition: vssolution3d.cpp:3428
gl3::GlBuilder::glVertex3d
void glVertex3d(double x, double y, double z)
Definition: types.hpp:332
GetLineWidth
float GetLineWidth()
Definition: aux_vis.cpp:1609
Normalize
int Normalize(double v[])
Definition: geom_utils.hpp:38
VisualizationSceneSolution3d::FindNewBox
virtual void FindNewBox(bool prepare)
Definition: vssolution3d.cpp:962
cut_TriGeoms
thread_local Array< int > cut_TriGeoms
Definition: vssolution3d.cpp:34
VisualizationSceneSolution3d::DrawRefinedPyramidLevelSurf
void DrawRefinedPyramidLevelSurf(gl3::GlDrawable &target, const DenseMatrix &verts, const Vector &vals, const int *RG, const int np, const int face_splits, const DenseMatrix *grad=NULL)
Definition: vssolution3d.cpp:3291
VisualizationSceneSolution3d::LiftRefinedSurf
void LiftRefinedSurf(int n, DenseMatrix &pointmat, Vector &values, int *RG)
Definition: vssolution3d.cpp:1267
wnd
thread_local SdlWindow * wnd
Definition: aux_vis.cpp:50
VisualizationSceneSolution3d::GetPyramidFaceSplits
static int GetPyramidFaceSplits(const Array< bool > &quad_diag, const Array< int > &faces, const Array< int > &ofaces)
Definition: vssolution3d.cpp:3277
GetLineWidthMS
float GetLineWidthMS()
Definition: aux_vis.cpp:1614
SetLineWidthMS
void SetLineWidthMS(float width_ms)
Definition: aux_vis.cpp:1599
cut_QuadPts
thread_local IntegrationRule cut_QuadPts
Definition: vssolution3d.cpp:31
gl3::GlBuilder::glEnd
void glEnd()
Definition: types.hpp:309
VisualizationSceneSolution3d::ToggleAttributes
virtual void ToggleAttributes(Array< int > &attr_list)
Definition: vssolution3d.cpp:933
VisualizationSceneSolution3d::GetFaceNormals
void GetFaceNormals(const int FaceNo, const int side, const IntegrationRule &ir, DenseMatrix &normals)
Definition: vssolution3d.cpp:1157
VisualizationSceneScalarData::GetSceneObjs
virtual gl3::SceneInfo GetSceneObjs()
Definition: vsdata.cpp:984
gl3::GlDrawable::addLine
void addLine(const Vert &v1, const Vert &v2)
Definition: types.hpp:672
VisualizationSceneSolution3d::UpdateValueRange
virtual void UpdateValueRange(bool prepare)
Definition: vssolution3d.cpp:1056
gl3::Vertex::create
static Vertex create(double *d)
Definition: types.hpp:207
VisualizationSceneScalarData::ComputeBdrAttrCenter
void ComputeBdrAttrCenter()
Compute the center of gravity for each boundary attribute.
Definition: vsdata.cpp:1612
gl3::GlBuilder::glNormal3dv
void glNormal3dv(const double *d)
Definition: types.hpp:407
Plane
Definition: vsdata.hpp:25
glTF_Builder::addBuffer
buffer_id addBuffer(const std::string &bufferName)
Definition: gltf.cpp:24
VisualizationSceneSolution3d::DrawRefinedSurfLevelLines
void DrawRefinedSurfLevelLines(int n, DenseMatrix &pointmat, Vector &values, Array< int > &RefGeoms)
Definition: vssolution3d.cpp:1397
gl3::GlBuilder::glBegin
void glBegin(GLenum e)
Definition: types.hpp:295
Compute3DUnitNormal
int Compute3DUnitNormal(const double p1[], const double p2[], const double p3[], double nor[])
Definition: geom_utils.hpp:101
VisualizationSceneSolution3d::PrepareOrderingCurve1
virtual void PrepareOrderingCurve1(gl3::GlDrawable &buf, bool arrows, bool color)
Definition: vssolution3d.cpp:143
VisualizationSceneSolution3d::DrawRefinedWedgeLevelSurf
void DrawRefinedWedgeLevelSurf(gl3::GlDrawable &target, const DenseMatrix &verts, const Vector &vals, const int *RG, const int np, const int face_splits, const DenseMatrix *grad=NULL)
Definition: vssolution3d.cpp:3445
VisualizationSceneSolution3d::DrawRefinedSurfEdges
void DrawRefinedSurfEdges(int n, DenseMatrix &pointmat, Vector &values, Array< int > &RefEdges, int part=-1)
Definition: vssolution3d.cpp:1368
CutReferenceElements
void CutReferenceElements(int TimesToRefine, double lambda)
Definition: vssolution3d.cpp:1609
GetMultisample
int GetMultisample()
Definition: aux_vis.cpp:1573
glTF_Builder
Definition: gltf.hpp:24
vssol3d
thread_local VisualizationSceneSolution3d * vssol3d
Definition: vssolution3d.cpp:24
VisualizationSceneScalarData::shrinkmat
double shrinkmat
Shrink factor with respect to the element (material) attributes centers.
Definition: vsdata.hpp:136
vs
thread_local VisualizationSceneScalarData * vs
Definition: glvis.cpp:67
gl3::GlBuilder::glNormal3d
void glNormal3d(double nx, double ny, double nz)
Definition: types.hpp:401
gl3::GlDrawable::createBuilder
GlBuilder createBuilder()
Definition: types.hpp:731
cut_QuadGeoms
thread_local Array< int > cut_QuadGeoms
Definition: vssolution3d.cpp:32
VisualizationSceneSolution3d::GetHexFaceSplits
static int GetHexFaceSplits(const Array< bool > &quad_diag, const Array< int > &faces, const Array< int > &ofaces)
Definition: vssolution3d.cpp:3574
glTF_Builder::addMesh
mesh_id addMesh(const std::string &meshName)
Definition: gltf.cpp:322
InnerProd
double InnerProd(const double a[], const double b[])
Definition: geom_utils.hpp:26
VisualizationSceneSolution3d::DrawTetLevelSurf
void DrawTetLevelSurf(gl3::GlDrawable &target, const DenseMatrix &verts, const Vector &vals, const int *ind, const Array< double > &levels, const DenseMatrix *grad=NULL)
Definition: vssolution3d.cpp:3109
VisualizationSceneSolution3d::SetShading
virtual void SetShading(int, bool)
Definition: vssolution3d.cpp:847
VisualizationSceneSolution3d::DrawRefinedSurf
void DrawRefinedSurf(int n, double *points, int elem, int func, int part=-1)
Definition: vssolution3d.cpp:1200
SetLineWidth
void SetLineWidth(float width)
Definition: aux_vis.cpp:1590
triangle_counter
thread_local int triangle_counter
Definition: vssolution3d.cpp:3106
VisualizationSceneSolution3d::GetSceneObjs
virtual gl3::SceneInfo GetSceneObjs()
Definition: vssolution3d.cpp:3940
VisualizationSceneSolution3d::CutRefinedFace
void CutRefinedFace(gl3::GlDrawable &target, const DenseMatrix &verts, const Vector &vert_dist, const Vector &vals, const Geometry::Type geom, const int *faces, int num_faces)
Definition: vssolution3d.cpp:2745
gl3::GlBuilder
Crude fixed-function OpenGL emulation helper.
Definition: types.hpp:261
glTF_Builder::writeFile
int writeFile()
Definition: gltf.cpp:472
KeyiPressed
void KeyiPressed()
Definition: vssolution.cpp:317
glTF_Builder::addNodeMesh
void addNodeMesh(node_id node, mesh_id mesh)
Definition: gltf.cpp:433
VisualizationSceneSolution3d::CutRefinedElement
void CutRefinedElement(gl3::GlDrawable &target, const DenseMatrix &verts, const Vector &vert_dist, const Vector &vals, const Geometry::Type geom, const int *elems, int num_elems, int func)
Definition: vssolution3d.cpp:2618
KeyIPressed
void KeyIPressed()
Definition: vssolution.cpp:323
ToggleMagicKey
void ToggleMagicKey()
Definition: vssolution3d.cpp:442
VisualizationSceneSolution3d::DrawRefinedHexLevelSurf
void DrawRefinedHexLevelSurf(gl3::GlDrawable &target, const DenseMatrix &verts, const Vector &vals, const int *RG, const int nh, const int face_splits, const DenseMatrix *grad=NULL)
Definition: vssolution3d.cpp:3591
VisualizationSceneSolution3d::PrepareOrderingCurve
virtual void PrepareOrderingCurve()
Definition: vssolution3d.cpp:131
KeyUPressed
void KeyUPressed()
Definition: vsvector.cpp:203
cut_TriPts
thread_local IntegrationRule cut_TriPts
Definition: vssolution3d.cpp:33
VisualizationSceneSolution3d::SetRefineFactors
virtual void SetRefineFactors(int, int)
Definition: vssolution3d.cpp:889
gl3::GlBuilder::glVertex3dv
void glVertex3dv(const double *d)
Definition: types.hpp:396
gl3::RenderParams::contains_translucent
bool contains_translucent
Definition: renderer.hpp:55
KeyVPressed
void KeyVPressed()
Definition: vsvector.cpp:147
KeyuPressed
void KeyuPressed()
Definition: vsvector.cpp:158
gl3::RenderParams::static_color
std::array< float, 4 > static_color
Definition: renderer.hpp:48
VisualizationSceneSolution3d::GetHelpString
virtual std::string GetHelpString() const
Definition: vssolution3d.cpp:38
SdlWindow::setOnKeyDown
void setOnKeyDown(int key, Delegate func)
Definition: sdl.hpp:188
gl3::SceneInfo::queue
RenderQueue queue
Definition: renderer.hpp:63
VisualizationSceneScalarData::ComputeElemAttrCenter
void ComputeElemAttrCenter()
Compute the center of gravity for each element attribute.
Definition: vsdata.cpp:1641
VisualizationSceneSolution3d::FindNewValueRange
virtual void FindNewValueRange(bool prepare)
Definition: vssolution3d.cpp:1028
gl3::RenderParams::mesh_material
Material mesh_material
Definition: renderer.hpp:44
gl3::RenderParams::clip_plane_eqn
std::array< double, 4 > clip_plane_eqn
Definition: renderer.hpp:52
VisualizationSceneScalarData::DoAutoscale
void DoAutoscale(bool prepare)
Definition: vsdata.cpp:48
VisualizationSceneSolution3d::NewMeshAndSolution
void NewMeshAndSolution(Mesh *new_m, Vector *new_sol, GridFunction *new_u=NULL)
Definition: vssolution3d.cpp:812
VisualizationScene::BLK_MAT
const Material BLK_MAT
Definition: openglvis.hpp:78
Generated on Tue May 24 2022 09:49:54 for GLVis by   doxygen 1.8.5