GLVis  v4.2
Accurate and flexible finite element visualization
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
vsvector3d.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 
20 using namespace mfem;
21 using namespace std;
22 
23 // Reference geometry with a cut in the middle, which subdivides GeometryRefiner
24 // when cut_lambda is updated, see keys Ctrl+F3/F4. These variables are defined
25 // in lib/vssolution3d.cpp.
26 extern thread_local IntegrationRule cut_QuadPts;
27 extern thread_local Array<int> cut_QuadGeoms;
28 extern thread_local IntegrationRule cut_TriPts;
29 extern thread_local Array<int> cut_TriGeoms;
30 extern void CutReferenceElements(int TimesToRefine, double lambda);
31 
32 
33 std::string VisualizationSceneVector3d::GetHelpString() const
34 {
35  std::stringstream os;
36  os << endl
37  << "+------------------------------------+" << endl
38  << "| Keys |" << endl
39  << "+------------------------------------+" << endl
40  << "| a - Displays/Hides the axes |" << endl
41  << "| A - Turns antialiasing on/off |" << endl
42  << "| b - Displacements step back |" << endl
43  << "| c - Toggle colorbar and caption |" << endl
44  << "| C - Change the main plot caption |" << endl
45  << "| d - Displays/Hides displacements |" << endl
46  << "| e - Displays/Hides the elements |" << endl
47  << "| E - Toggle the elements in the CP |" << endl
48  << "| f - Smooth/Flat shading |" << endl
49  << "| F - Display mag./x/y/z component |" << endl
50  << "| g - Toggle background |" << endl
51  << "| h - Displays help menu |" << endl
52  << "| i - Toggle cutting plane |" << endl
53  << "| j - Turn on/off perspective |" << endl
54  << "| k/K Adjust the transparency level |" << endl
55  << "| ,/< Adjust color transparency |" << endl
56  << "| l - Turns on/off the light |" << endl
57  << "| L - Toggle logarithmic scale |" << endl
58  << "| m - Displays/Hides the mesh |" << endl
59  << "| M - Toggle the mesh in the CP |" << endl
60  << "| n - Displacements step forward |" << endl
61  << "| o/O (De)refine elem, disc shading |" << endl
62  << "| p/P Cycle through color palettes |" << endl
63  << "| q - Quits |" << endl
64  << "| r - Reset the plot to 3D view |" << endl
65  << "| R - Reset the plot to 2D view |" << endl
66  << "| s - Turn on/off unit cube scaling |" << endl
67  << "| S - Take snapshot/Record a movie |" << endl
68  << "| t - Cycle materials and lights |" << endl
69  << "| u/U Move the level field vectors |" << endl
70  << "| v/V Vector field |" << endl
71  << "| w/W Add/Delete level field vector |" << endl
72  << "| x/X Rotate cutting plane (phi) |" << endl
73  << "| y/Y Rotate cutting plane (theta) |" << endl
74  << "| z/Z Translate cutting plane |" << endl
75  << "| \\ - Set light source position |" << endl
76  << "| Ctrl+p - Print to a PDF file |" << endl
77  << "+------------------------------------+" << endl
78  << "| Function keys |" << endl
79  << "+------------------------------------+" << endl
80  << "| F1 - X window info and keystrokes |" << endl
81  << "| F2 - Update colors, etc. |" << endl
82  << "| F3/F4 - Shrink/Zoom bdr elements |" << endl
83  << "| Ctrl+F3/F4 - Cut face bdr elements |" << endl
84  << "| F6 - Palette options |" << endl
85  << "| F7 - Manually set min/max value |" << endl
86  << "| F8 - List of subdomains to show |" << endl
87  << "| F9/F10 - Walk through subdomains |" << endl
88  << "| F11/F12 - Shrink/Zoom subdomains |" << endl
89  << "+------------------------------------+" << endl
90  << "| Keypad |" << endl
91  << "+------------------------------------+" << endl
92  << "| 1-9 Small rotation, reset with 5 |" << endl
93  << "| *,/ Scale up/down |" << endl
94  << "| +/- Change z-scaling |" << endl
95  << "| . - Start/stop spinning |" << endl
96  << "| 0/Enter - Spinning speed and dir. |" << endl
97  << "+------------------------------------+" << endl
98  << "| Mouse |" << endl
99  << "+------------------------------------+" << endl
100  << "| left btn - Rotation |" << endl
101  << "| middle btn - Translation |" << endl
102  << "| right btn - Scaling |" << endl
103  << "| left + Alt - Tilt |" << endl
104  << "| left + Shift - Spinning |" << endl
105  << "| right + Shift - Change light pos. |" << endl
106  << "| left + Ctrl - Spherical rotation |" << endl
107  << "| middle+ Ctrl - Object translation |" << endl
108  << "| right + Ctrl - Object scaling |" << endl
109  << "| left + Ctrl + Shift - z-Spinning |" << endl
110  << "+------------------------------------+" << endl;
111  return os.str();
112 }
113 
114 thread_local VisualizationSceneVector3d *vsvector3d;
115 extern thread_local VisualizationScene *locscene;
116 extern thread_local GeometryRefiner GLVisGeometryRefiner;
117 
118 static void KeyDPressed()
119 {
120  vsvector3d -> ToggleDisplacements();
121  SendExposeEvent();
122 }
123 
124 static void KeyNPressed()
125 {
126  if (vsvector3d -> drawdisp)
127  vsvector3d -> ianimd =( (vsvector3d ->ianimd + 1) %
128  (vsvector3d -> ianimmax + 1) );
129  else
130  vsvector3d -> ianim =( (vsvector3d ->ianim + 1) %
131  (vsvector3d -> ianimmax + 1) );
132  vsvector3d -> NPressed();
133 }
134 
135 static void KeyBPressed()
136 {
137  if (vsvector3d -> drawdisp)
138  vsvector3d ->ianimd = ((vsvector3d ->ianimd +
139  vsvector3d ->ianimmax) %
140  (vsvector3d ->ianimmax + 1));
141  else
142  vsvector3d ->ianim = ((vsvector3d ->ianim +
143  vsvector3d ->ianimmax) %
144  (vsvector3d ->ianimmax + 1));
145  vsvector3d -> NPressed();
146 }
147 
148 static void KeyrPressed()
149 {
150  locscene -> spinning = 0;
151  RemoveIdleFunc(MainLoop);
152  vsvector3d -> CenterObject();
153  locscene -> ViewAngle = 45.0;
154  locscene -> ViewScale = 1.0;
155  locscene -> ViewCenterX = 0.0;
156  locscene -> ViewCenterY = 0.0;
157  vsvector3d -> ianim = vsvector3d -> ianimd = 0;
158  vsvector3d -> Prepare();
159  vsvector3d -> PrepareLines();
160  vsvector3d -> PrepareDisplacedMesh();
161  vsvector3d -> key_r_state = 0;
162  SendExposeEvent();
163 }
164 
165 static void KeyRPressed()
166 {
167  locscene->spinning = 0;
168  RemoveIdleFunc(MainLoop);
169  vsvector3d -> ianim = vsvector3d -> ianimd = 0;
170  vsvector3d -> Prepare();
171  vsvector3d -> PrepareLines();
172  vsvector3d -> PrepareDisplacedMesh();
173  vsvector3d->Toggle2DView();
174  SendExposeEvent();
175 }
176 
177 void VisualizationSceneVector3d::NPressed()
178 {
179  if (drawdisp)
180  {
181  PrepareDisplacedMesh();
182  }
183  else
184  {
185  Prepare();
186  PrepareLines();
187  }
188 
189  SendExposeEvent();
190 }
191 
192 static void KeyuPressed()
193 {
194  vsvector3d -> ToggleVectorFieldLevel(+1);
195  SendExposeEvent();
196 }
197 
198 static void KeyUPressed()
199 {
200  vsvector3d -> ToggleVectorFieldLevel(-1);
201  SendExposeEvent();
202 }
203 
204 void VisualizationSceneVector3d::ToggleVectorFieldLevel(int v)
205 {
206  int i;
207  for (i = 0; i < vflevel.Size(); i++)
208  if (vflevel[i] == 0 && v == -1)
209  {
210  vflevel[i] = nl;
211  }
212  else
213  {
214  vflevel[i] = (vflevel[i] + v) % (nl+1);
215  }
216  for (i = 0; i < vflevel.Size(); i++)
217  {
218  dvflevel[i] = level[vflevel[i]];
219  }
220  vsvector3d -> PrepareVectorField();
221 }
222 
223 static void KeywPressed()
224 {
225  vsvector3d -> AddVectorFieldLevel();
226  SendExposeEvent();
227 }
228 
229 static void KeyWPressed()
230 {
231  vsvector3d -> RemoveVectorFieldLevel();
232  SendExposeEvent();
233 }
234 
235 void VisualizationSceneVector3d::AddVectorFieldLevel()
236 {
237  int next = vflevel[vflevel.Size()-1];
238  next = (next + 1) % (nl+1);
239  vflevel.Append(next);
240  dvflevel.Append(level[next]);
241  vsvector3d -> PrepareVectorField();
242 }
243 
244 void VisualizationSceneVector3d::RemoveVectorFieldLevel()
245 {
246  vflevel.DeleteLast();
247  dvflevel.DeleteLast();
248  vsvector3d -> PrepareVectorField();
249 }
250 
251 static void KeyvPressed()
252 {
253  vsvector3d -> ToggleVectorField(1);
254  SendExposeEvent();
255 }
256 
257 static void KeyVPressed()
258 {
259  vsvector3d -> ToggleVectorField(-1);
260  SendExposeEvent();
261 }
262 
263 static void VectorKeyFPressed()
264 {
265  vsvector3d->ToggleScalarFunction();
266  SendExposeEvent();
267 }
268 
269 void VisualizationSceneVector3d::ToggleVectorField(int i)
270 {
271  drawvector = (drawvector+i+6)%6;
272  PrepareVectorField();
273 }
274 
275 static const char *scal_func_name[] =
276 {"magnitude", "x-component", "y-component", "z-component"};
277 
278 void VisualizationSceneVector3d::SetScalarFunction()
279 {
280  FiniteElementSpace *fes = (VecGridF) ? VecGridF->FESpace() : NULL;
281 
282  switch (scal_func)
283  {
284  case 0: // magnitude
285  for (int i = 0; i < sol->Size(); i++)
286  (*sol)(i) = sqrt((*solx)(i) * (*solx)(i) +
287  (*soly)(i) * (*soly)(i) +
288  (*solz)(i) * (*solz)(i) );
289  if (GridF)
290  {
291  Array<int> dofs(3);
292  for (int i = 0; i < GridF->Size(); i++)
293  {
294  dofs.SetSize(1);
295  dofs[0] = i;
296  fes->DofsToVDofs(dofs);
297  double x = (*VecGridF)(dofs[0]);
298  double y = (*VecGridF)(dofs[1]);
299  double z = (*VecGridF)(dofs[2]);
300 
301  (*GridF)(i) = sqrt(x*x+y*y+z*z);
302  }
303  }
304  break;
305  case 1: // x-component
306  *sol = *solx;
307  if (GridF)
308  for (int i = 0; i < GridF->Size(); i++)
309  {
310  (*GridF)(i) = (*VecGridF)(fes->DofToVDof(i, 0));
311  }
312  break;
313  case 2: // y-component
314  *sol = *soly;
315  if (GridF)
316  for (int i = 0; i < GridF->Size(); i++)
317  {
318  (*GridF)(i) = (*VecGridF)(fes->DofToVDof(i, 1));
319  }
320  break;
321  case 3: // z-component
322  *sol = *solz;
323  if (GridF)
324  for (int i = 0; i < GridF->Size(); i++)
325  {
326  (*GridF)(i) = (*VecGridF)(fes->DofToVDof(i, 2));
327  }
328  break;
329  }
330  extra_caption = scal_func_name[scal_func];
331 }
332 
333 void VisualizationSceneVector3d::ToggleScalarFunction()
334 {
335  scal_func = (scal_func + 1) % 4;
336  cout << "Displaying " << scal_func_name[scal_func] << endl;
337  SetScalarFunction();
338  FindNewValueRange(true);
339 }
340 
341 VisualizationSceneVector3d::VisualizationSceneVector3d(Mesh &m, Vector &sx,
342  Vector &sy, Vector &sz)
343 {
344  mesh = &m;
345  solx = &sx;
346  soly = &sy;
347  solz = &sz;
348 
349  sol = new Vector(mesh->GetNV());
350 
351  sfes = NULL;
352  VecGridF = NULL;
353 
354  Init();
355 }
356 
357 VisualizationSceneVector3d::VisualizationSceneVector3d(GridFunction &vgf)
358 {
359  FiniteElementSpace *fes = vgf.FESpace();
360  if (fes == NULL || fes->GetVDim() != 3)
361  {
362  cout << "VisualizationSceneVector3d::VisualizationSceneVector3d" << endl;
363  exit(1);
364  }
365 
366  VecGridF = &vgf;
367 
368  mesh = fes->GetMesh();
369 
370  sfes = new FiniteElementSpace(mesh, fes->FEColl(), 1, fes->GetOrdering());
371  GridF = new GridFunction(sfes);
372 
373  solx = new Vector(mesh->GetNV());
374  soly = new Vector(mesh->GetNV());
375  solz = new Vector(mesh->GetNV());
376 
377  vgf.GetNodalValues(*solx, 1);
378  vgf.GetNodalValues(*soly, 2);
379  vgf.GetNodalValues(*solz, 3);
380 
381  sol = new Vector(mesh->GetNV());
382 
383  Init();
384 }
385 
386 void VisualizationSceneVector3d::Init()
387 {
388  key_r_state = 0;
389 
390  drawdisp = 0;
391  drawvector = 0;
392  scal_func = 0;
393 
394  ianim = ianimd = 0;
395  ianimmax = 10;
396 
397  SetScalarFunction();
398 
399  VisualizationSceneSolution3d::Init();
400 
401  PrepareVectorField();
402  PrepareDisplacedMesh();
403 
404  vflevel.Append(0);
405  dvflevel.Append(level[0]);
406 
407  vsvector3d = this;
408 
409  // static int init = 0;
410  // if (!init)
411  {
412  // init = 1;
413 
414  wnd->setOnKeyDown('d', KeyDPressed);
415  wnd->setOnKeyDown('D', KeyDPressed);
416 
417  wnd->setOnKeyDown('n', KeyNPressed);
418  wnd->setOnKeyDown('N', KeyNPressed);
419 
420  wnd->setOnKeyDown('b', KeyBPressed);
421  wnd->setOnKeyDown('B', KeyBPressed);
422 
423  wnd->setOnKeyDown('r', KeyrPressed); // adds another function to 'r' and 'R'
424  wnd->setOnKeyDown('R', KeyRPressed); // the other function is in vsdata.cpp
425 
426  wnd->setOnKeyDown('u', KeyuPressed); // Keys u, U are also used in
427  wnd->setOnKeyDown('U', KeyUPressed); // VisualizationSceneSolution3d
428 
429  wnd->setOnKeyDown('w', KeywPressed); // Keys w, W are also used in
430  wnd->setOnKeyDown('W', KeyWPressed); // VisualizationSceneSolution3d
431 
432  wnd->setOnKeyDown('v', KeyvPressed); // Keys v, V are also used in
433  wnd->setOnKeyDown('V', KeyVPressed); // VisualizationSceneSolution3d
434 
435  wnd->setOnKeyDown('F', VectorKeyFPressed);
436  }
437 }
438 
439 VisualizationSceneVector3d::~VisualizationSceneVector3d()
440 {
441  delete sol;
442 
443  if (VecGridF)
444  {
445  delete solz;
446  delete soly;
447  delete solx;
448  delete GridF;
449  delete sfes;
450  }
451 }
452 
453 void VisualizationSceneVector3d::NewMeshAndSolution(
454  Mesh *new_m, GridFunction *new_v)
455 {
456  delete sol;
457  if (VecGridF)
458  {
459  delete solz;
460  delete soly;
461  delete solx;
462  delete GridF;
463  delete sfes;
464  }
465  if (mesh->GetNV() != new_m->GetNV())
466  {
467  delete [] node_pos;
468  node_pos = new double[new_m->GetNV()];
469  }
470 
471  // If the number of surface elements changes, recompute the refinement factor
472  if (mesh->Dimension() != new_m->Dimension() ||
473  (mesh->Dimension() == 2 && mesh->GetNE() != new_m->GetNE()) ||
474  (mesh->Dimension() == 3 && mesh->GetNBE() != new_m->GetNBE()))
475  {
476  mesh = new_m;
477  int ref = GetAutoRefineFactor();
478  if (TimesToRefine != ref)
479  {
480  TimesToRefine = ref;
481  cout << "Subdivision factor = " << TimesToRefine << endl;
482  }
483  }
484 
485  FiniteElementSpace *new_fes = new_v->FESpace();
486 
487  VecGridF = new_v;
488  mesh = new_m;
489  FindNodePos();
490 
491  sfes = new FiniteElementSpace(mesh, new_fes->FEColl(), 1,
492  new_fes->GetOrdering());
493  GridF = new GridFunction(sfes);
494 
495  solx = new Vector(mesh->GetNV());
496  soly = new Vector(mesh->GetNV());
497  solz = new Vector(mesh->GetNV());
498 
499  VecGridF->GetNodalValues(*solx, 1);
500  VecGridF->GetNodalValues(*soly, 2);
501  VecGridF->GetNodalValues(*solz, 3);
502 
503  sol = new Vector(mesh->GetNV());
504 
505  SetScalarFunction();
506 
507  DoAutoscale(false);
508 
509  Prepare();
510  PrepareLines();
511  CPPrepare();
512  PrepareLevelSurf();
513 
514  PrepareVectorField();
515  PrepareDisplacedMesh();
516 }
517 
518 void VisualizationSceneVector3d::PrepareFlat()
519 {
520  int i, j;
521 
522  disp_buf.clear();
523 
524  int dim = mesh->Dimension();
525  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
526  DenseMatrix pointmat;
527  Array<int> vertices;
528  double p[4][3], c[4];
529 
530  for (i = 0; i < ne; i++)
531  {
532  if (dim == 3)
533  {
534  if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) { continue; }
535 
536  mesh->GetBdrPointMatrix(i, pointmat);
537  mesh->GetBdrElementVertices(i, vertices);
538  }
539  else
540  {
541  if (!bdr_attr_to_show[mesh->GetAttribute(i)-1]) { continue; }
542 
543  mesh->GetPointMatrix(i, pointmat);
544  mesh->GetElementVertices(i, vertices);
545  }
546 
547  for (j = 0; j < pointmat.Width(); j++)
548  {
549  pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax;
550  pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax;
551  pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax;
552  }
553 
554  for (j = 0; j < pointmat.Width(); j++)
555  {
556  p[j][0] = pointmat(0, j);
557  p[j][1] = pointmat(1, j);
558  p[j][2] = pointmat(2, j);
559  c[j] = (*sol)(vertices[j]);
560  }
561  if (j == 3)
562  {
563  if (cut_lambda > 0)
564  {
565  DrawCutTriangle(disp_buf, p, c, minv, maxv);
566  }
567  else
568  {
569  DrawTriangle(disp_buf, p, c, minv, maxv);
570  }
571  }
572  else
573  {
574  if (cut_lambda > 0)
575  {
576  DrawCutQuad(disp_buf, p, c, minv, maxv);
577  }
578  else
579  {
580  DrawQuad(disp_buf, p, c, minv, maxv);
581  }
582  }
583  }
584  updated_bufs.emplace_back(&disp_buf);
585 }
586 
587 void VisualizationSceneVector3d::PrepareFlat2()
588 {
589  int fn, fo, di = 0, have_normals;
590  double bbox_diam, vmin, vmax;
591  int dim = mesh->Dimension();
592  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
593 
594  DenseMatrix pointmat, normals;
595  DenseMatrix vec_vals;
596  Vector values, normal;
597  RefinedGeometry * RefG;
598  Array<int> vertices;
599  double norm[3];
600  IsoparametricTransformation T;
601 
602  bbox_diam = sqrt( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
603  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
604  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
605  double sc = FaceShiftScale * bbox_diam;
606 
607  disp_buf.clear();
608 
609  vmin = numeric_limits<double>::infinity();
610  vmax = -vmin;
611  for (int i = 0; i < ne; i++)
612  {
613  int sides;
614  switch ((dim == 3) ? mesh->GetBdrElementType(i) : mesh->GetElementType(i))
615  {
616  case Element::TRIANGLE:
617  sides = 3;
618  break;
619 
620  case Element::QUADRILATERAL:
621  default:
622  sides = 4;
623  break;
624  }
625 
626  if (dim == 3)
627  {
628  if (!bdr_attr_to_show[mesh->GetBdrAttribute(i)-1]) { continue; }
629 
630  if (cplane == 2)
631  {
632  // for cplane == 2, get vertices of the volume element, not bdr
633  int f, o, e1, e2;
634  mesh->GetBdrElementFace(i, &f, &o);
635  mesh->GetFaceElements(f, &e1, &e2);
636  mesh->GetElementVertices(e1, vertices);
637  }
638  else
639  {
640  mesh->GetBdrElementVertices(i, vertices);
641  }
642  }
643  else
644  {
645  if (!bdr_attr_to_show[mesh->GetAttribute(i)-1]) { continue; }
646  mesh->GetElementVertices(i, vertices);
647  }
648 
649  if (cplane == 2 && CheckPositions(vertices)) { continue; }
650 
651  if (dim == 3)
652  {
653  mesh -> GetBdrElementFace (i, &fn, &fo);
654  RefG = GLVisGeometryRefiner.Refine(mesh -> GetFaceBaseGeometry (fn),
655  TimesToRefine);
656  if (!cut_updated)
657  {
658  // Update the cut version of the reference geometries
659  CutReferenceElements(TimesToRefine, cut_lambda);
660  cut_updated = true;
661  }
662 
663  // di = GridF->GetFaceValues(fn, 2, RefG->RefPts, values, pointmat);
664  di = fo % 2;
665  if (di == 1 && !mesh->FaceIsInterior(fn))
666  {
667  di = 0;
668  }
669 
670  IntegrationRule &RefPts = (cut_lambda > 0) ?
671  ((sides == 3) ? cut_TriPts : cut_QuadPts) :
672  RefG->RefPts;
673  GridF->GetFaceValues(fn, di, RefPts, values, pointmat);
674  if (ianim > 0)
675  {
676  VecGridF->GetFaceVectorValues(fn, di, RefPts, vec_vals,
677  pointmat);
678  pointmat.Add(double(ianim)/ianimmax, vec_vals);
679  have_normals = 0;
680  }
681  else
682  {
683  GetFaceNormals(fn, di, RefPts, normals);
684  have_normals = 1;
685  }
686  ShrinkPoints(pointmat, i, fn, di);
687  }
688  else // dim == 2
689  {
690  RefG = GLVisGeometryRefiner.Refine(mesh->GetElementBaseGeometry(i),
691  TimesToRefine);
692  if (!cut_updated)
693  {
694  // Update the cut version of the reference geometries
695  CutReferenceElements(TimesToRefine, cut_lambda);
696  cut_updated = true;
697  }
698  IntegrationRule &RefPts = (cut_lambda > 0) ?
699  ((sides == 3) ? cut_TriPts : cut_QuadPts) :
700  RefG->RefPts;
701  GridF->GetValues(i, RefPts, values, pointmat);
702  if (ianim > 0)
703  {
704  VecGridF->GetVectorValues(i, RefPts, vec_vals, pointmat);
705  pointmat.Add(double(ianim)/ianimmax, vec_vals);
706  have_normals = 0;
707  }
708  else
709  {
710  const IntegrationRule &ir = (cut_lambda > 0) ?
711  ((sides == 3) ? cut_TriPts : cut_QuadPts) :
712  RefG->RefPts;
713  normals.SetSize(3, values.Size());
714  mesh->GetElementTransformation(i, &T);
715  for (int j = 0; j < values.Size(); j++)
716  {
717  T.SetIntPoint(&ir.IntPoint(j));
718  normals.GetColumnReference(j, normal);
719  CalcOrtho(T.Jacobian(), normal);
720  normal /= normal.Norml2();
721  }
722  have_normals = 1;
723  di = 0;
724  }
725  ShrinkPoints(pointmat, i, 0, 0);
726  }
727 
728  vmin = fmin(vmin, values.Min());
729  vmax = fmax(vmax, values.Max());
730 
731  if (sc != 0.0 && have_normals)
732  {
733  for (int j = 0; j < 3; j++)
734  {
735  norm[j] = 0.0;
736  }
737  Normalize(normals);
738  for (int k = 0; k < normals.Width(); k++)
739  for (int j = 0; j < 3; j++)
740  {
741  norm[j] += normals(j, k);
742  }
743  Normalize(norm);
744  for (int k = 0; k < pointmat.Width(); k++)
745  {
746  double val = sc * (values(k) - minv) / (maxv - minv);
747  for (int j = 0; j < 3; j++)
748  {
749  pointmat(j, k) += val * norm[j];
750  }
751  }
752  have_normals = 0;
753  }
754 
755  have_normals = have_normals ? 2 : 0;
756  if (di)
757  {
758  have_normals = -1 - have_normals;
759  }
760 
761  Array<int> &RefGeoms = (cut_lambda > 0) ?
762  ((sides == 3) ? cut_TriGeoms : cut_QuadGeoms) :
763  RefG->RefGeoms;
764  int psides = (cut_lambda > 0) ? 4 : sides;
765  DrawPatch(disp_buf, pointmat, values, normals, psides, RefGeoms,
766  minv, maxv, have_normals);
767  }
768  updated_bufs.emplace_back(&disp_buf);
769  cout << "VisualizationSceneVector3d::PrepareFlat2() : [min,max] = ["
770  << vmin << "," << vmax << "]" << endl;
771 }
772 
773 void VisualizationSceneVector3d::Prepare()
774 {
775  int i,j;
776 
777  switch (shading)
778  {
779  case 0:
780  PrepareFlat();
781  return;
782  case 2:
783  PrepareFlat2();
784  return;
785  default:
786  break;
787  }
788 
789  disp_buf.clear();
790  gl3::GlBuilder draw = disp_buf.createBuilder();
791 
792  int dim = mesh->Dimension();
793  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
794  int nv = mesh -> GetNV();
795  DenseMatrix pointmat;
796  Array<int> vertices;
797  double nor[3];
798 
799  Vector nx(nv);
800  Vector ny(nv);
801  Vector nz(nv);
802 
803  const Array<int> &attributes =
804  ((dim == 3) ? mesh->bdr_attributes : mesh->attributes);
805  for (int d = 0; d < attributes.Size(); d++)
806  {
807  if (!bdr_attr_to_show[attributes[d]-1]) { continue; }
808 
809  nx=0.;
810  ny=0.;
811  nz=0.;
812 
813  for (i = 0; i < ne; i++)
814  {
815  int attr =
816  (dim == 3) ? mesh->GetBdrAttribute(i) : mesh->GetAttribute(i);
817  if (attr != attributes[d]) { continue; }
818 
819  if (dim == 3)
820  {
821  mesh->GetBdrPointMatrix (i, pointmat);
822  mesh->GetBdrElementVertices (i, vertices);
823  }
824  else
825  {
826  mesh->GetPointMatrix(i, pointmat);
827  mesh->GetElementVertices(i, vertices);
828  }
829 
830  for (j = 0; j < pointmat.Size(); j++)
831  {
832  pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax;
833  pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax;
834  pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax;
835  }
836 
837  if (pointmat.Width() == 3)
838  j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1),
839  &pointmat(0,2), nor);
840  else
841  j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1),
842  &pointmat(0,2), &pointmat(0,3), nor);
843  if (j == 0)
844  for (j = 0; j < pointmat.Size(); j++)
845  {
846  nx(vertices[j]) += nor[0];
847  ny(vertices[j]) += nor[1];
848  nz(vertices[j]) += nor[2];
849  }
850  }
851 
852  for (i = 0; i < ne; i++)
853  {
854  int attr =
855  (dim == 3) ? mesh->GetBdrAttribute(i) : mesh->GetAttribute(i);
856  if (attr != attributes[d]) { continue; }
857 
858  int el_type =
859  (dim == 3) ? mesh->GetBdrElementType(i) : mesh->GetElementType(i);
860  switch (el_type)
861  {
862  case Element::TRIANGLE:
863  draw.glBegin (GL_TRIANGLES);
864  break;
865 
866  case Element::QUADRILATERAL:
867  draw.glBegin (GL_QUADS);
868  break;
869  }
870  if (dim == 3)
871  {
872  mesh->GetBdrPointMatrix (i, pointmat);
873  mesh->GetBdrElementVertices (i, vertices);
874  }
875  else
876  {
877  mesh->GetPointMatrix(i, pointmat);
878  mesh->GetElementVertices(i, vertices);
879  }
880  for (j = 0; j < pointmat.Size(); j++)
881  {
882  pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax;
883  pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax;
884  pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax;
885  }
886  for (j = 0; j < pointmat.Size(); j++)
887  {
888  MySetColor(draw, (*sol)(vertices[j]), minv, maxv);
889  draw.glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j]));
890  draw.glVertex3dv(&pointmat(0, j));
891  }
892  draw.glEnd();
893  }
894  }
895  updated_bufs.emplace_back(&disp_buf);
896 }
897 
898 void VisualizationSceneVector3d::PrepareLines()
899 {
900  if (!drawmesh) { return; }
901 
902  if (shading == 2)
903  {
904  PrepareLines2();
905  return;
906  }
907 
908  int dim = mesh->Dimension();
909  int i, j, ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
910  DenseMatrix pointmat;
911  Array<int> vertices;
912  double point[4][4];
913 
914  line_buf.clear();
915 
916  for (i = 0; i < ne; i++)
917  {
918  int attr = (dim == 3) ? mesh->GetBdrAttribute(i) : mesh->GetAttribute(i);
919  if (!bdr_attr_to_show[attr-1]) { continue; }
920 
921  if (dim == 3)
922  {
923  if (cplane == 2)
924  {
925  // for cplane == 2, get vertices of the volume element, not bdr
926  int f, o, e1, e2;
927  mesh->GetBdrElementFace(i, &f, &o);
928  mesh->GetFaceElements(f, &e1, &e2);
929  mesh->GetElementVertices(e1, vertices);
930 
931  if (CheckPositions(vertices)) { continue; }
932  }
933  mesh->GetBdrElementVertices(i, vertices);
934  mesh->GetBdrPointMatrix(i, pointmat);
935  }
936  else
937  {
938  mesh->GetElementVertices(i, vertices);
939  mesh->GetPointMatrix(i, pointmat);
940  }
941 
942  if (cplane == 2 && CheckPositions(vertices)) { continue; }
943 
944  for (j = 0; j < pointmat.Size(); j++)
945  {
946  pointmat(0, j) += (*solx)(vertices[j])*(ianim)/ianimmax;
947  pointmat(1, j) += (*soly)(vertices[j])*(ianim)/ianimmax;
948  pointmat(2, j) += (*solz)(vertices[j])*(ianim)/ianimmax;
949  }
950  gl3::GlBuilder line = line_buf.createBuilder();
951  switch (drawmesh)
952  {
953  case 1:
954  line.glBegin(GL_LINE_LOOP);
955 
956  for (j = 0; j < pointmat.Size(); j++)
957  {
958  line.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j));
959  }
960  line.glEnd();
961  break;
962 
963  case 2:
964  for (j = 0; j < pointmat.Size(); j++)
965  {
966  for (int k = 0; k < 3; k++)
967  {
968  point[j][k] = pointmat(k,j);
969  }
970  point[j][3] = (*sol)(vertices[j]);
971  }
972  DrawPolygonLevelLines(line, point[0], pointmat.Size(), level, false);
973  break;
974  }
975  }
976  updated_bufs.emplace_back(&line_buf);
977 }
978 
979 void VisualizationSceneVector3d::PrepareLines2()
980 {
981  int i, j, k, fn, fo, di = 0;
982  double bbox_diam;
983 
984  line_buf.clear();
985  gl3::GlBuilder line = line_buf.createBuilder();
986 
987  int dim = mesh->Dimension();
988  int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
989  DenseMatrix pointmat;
990  DenseMatrix vec_vals;
991  Vector values;
992  RefinedGeometry * RefG;
993  Array<int> vertices;
994 
995  bbox_diam = sqrt ( (bb.x[1]-bb.x[0])*(bb.x[1]-bb.x[0]) +
996  (bb.y[1]-bb.y[0])*(bb.y[1]-bb.y[0]) +
997  (bb.z[1]-bb.z[0])*(bb.z[1]-bb.z[0]) );
998  double sc = FaceShiftScale * bbox_diam;
999 
1000  for (i = 0; i < ne; i++)
1001  {
1002  int attr = (dim == 3) ? mesh->GetBdrAttribute(i) : mesh->GetAttribute(i);
1003  if (!bdr_attr_to_show[attr-1]) { continue; }
1004 
1005  if (dim == 3)
1006  {
1007  if (cplane == 2)
1008  {
1009  // for cplane == 2, get vertices of the volume element, not bdr
1010  int f, o, e1, e2;
1011  mesh->GetBdrElementFace(i, &f, &o);
1012  mesh->GetFaceElements(f, &e1, &e2);
1013  mesh->GetElementVertices(e1, vertices);
1014  }
1015  else
1016  {
1017  mesh->GetBdrElementVertices (i, vertices);
1018  }
1019  }
1020  else
1021  {
1022  mesh->GetElementVertices(i, vertices);
1023  }
1024 
1025  if (cplane == 2 && CheckPositions(vertices)) { continue; }
1026 
1027  if (dim == 3)
1028  {
1029  mesh -> GetBdrElementFace (i, &fn, &fo);
1030  RefG = GLVisGeometryRefiner.Refine(mesh -> GetFaceBaseGeometry (fn),
1031  TimesToRefine);
1032  // di = GridF->GetFaceValues(fn, 2, RefG->RefPts, values, pointmat);
1033  di = fo % 2;
1034  if (di == 1 && !mesh->FaceIsInterior(fn))
1035  {
1036  di = 0;
1037  }
1038  GridF -> GetFaceValues (fn, di, RefG->RefPts, values, pointmat);
1039  VecGridF -> GetFaceVectorValues (fn, di, RefG->RefPts,
1040  vec_vals, pointmat);
1041  if (ianim > 0)
1042  {
1043  pointmat.Add (double(ianim)/ianimmax, vec_vals);
1044  }
1045  ShrinkPoints(pointmat, i, fn, di);
1046  }
1047  else
1048  {
1049  RefG = GLVisGeometryRefiner.Refine(mesh->GetElementBaseGeometry(i),
1050  TimesToRefine);
1051  GridF->GetValues(i, RefG->RefPts, values, pointmat);
1052  VecGridF->GetVectorValues(i, RefG->RefPts, vec_vals, pointmat);
1053  if (ianim > 0)
1054  {
1055  pointmat.Add(double(ianim)/ianimmax, vec_vals);
1056  }
1057  ShrinkPoints(pointmat, i, 0, 0);
1058  }
1059 
1060  int *RG = &(RefG->RefGeoms[0]);
1061  double pts[][3] =
1062  {
1063  { pointmat(0,RG[0]), pointmat(1,RG[0]), pointmat(2,RG[0]) },
1064  { pointmat(0,RG[1]), pointmat(1,RG[1]), pointmat(2,RG[1]) },
1065  { pointmat(0,RG[2]), pointmat(1,RG[2]), pointmat(2,RG[2]) }
1066  };
1067  double norm[3];
1068  Compute3DUnitNormal (pts[0], pts[1], pts[2], norm);
1069  if (di != 0 && sc != 0.0)
1070  {
1071  norm[0] = -norm[0];
1072  norm[1] = -norm[1];
1073  norm[2] = -norm[2];
1074  }
1075 
1076  if (drawmesh == 1)
1077  {
1078  Array<int> &REdges = RefG->RefEdges;
1079  line.glBegin (GL_LINES);
1080  for (k = 0; k < REdges.Size()/2; k++)
1081  {
1082  int *RE = &(REdges[2*k]);
1083 
1084  if (sc == 0.0)
1085  {
1086  for (j = 0; j < 2; j++)
1087  line.glVertex3d (pointmat(0, RE[j]), pointmat(1, RE[j]),
1088  pointmat(2, RE[j]));
1089  }
1090  else
1091  {
1092  for (j = 0; j < 2; j++)
1093  {
1094  double val = sc * (values(RE[j]) - minv) / (maxv - minv);
1095  line.glVertex3d (pointmat(0, RE[j])+val*norm[0],
1096  pointmat(1, RE[j])+val*norm[1],
1097  pointmat(2, RE[j])+val*norm[2]);
1098  }
1099  }
1100  }
1101  line.glEnd();
1102  }
1103  else if (drawmesh == 2)
1104  {
1105  double point[4][4];
1106  int sides;
1107  switch ((dim == 3) ? mesh->GetBdrElementType(i) :
1108  mesh->GetElementType(i))
1109  {
1110  case Element::TRIANGLE:
1111  sides = 3;
1112  break;
1113 
1114  case Element::QUADRILATERAL:
1115  default:
1116  sides = 4;
1117  break;
1118  }
1119  for (k = 0; k < RefG->RefGeoms.Size()/sides; k++)
1120  {
1121  RG = &(RefG->RefGeoms[k*sides]);
1122 
1123  if (sc == 0.0)
1124  {
1125  for (j = 0; j < sides; j++)
1126  {
1127  for (int ii = 0; ii < 3; ii++)
1128  {
1129  point[j][ii] = pointmat(ii, RG[j]);
1130  }
1131  point[j][3] = values(RG[j]);
1132  }
1133  }
1134  else
1135  {
1136  for (j = 0; j < sides; j++)
1137  {
1138  double val = (values(RG[j]) - minv) / (maxv - minv);
1139  val *= sc;
1140  for (int ii = 0; ii < 3; ii++)
1141  {
1142  point[j][ii] = pointmat(ii, RG[j])+val*norm[ii];
1143  }
1144  point[j][3] = values(RG[j]);
1145  }
1146  }
1147  DrawPolygonLevelLines(line, point[0], sides, level, false);
1148  }
1149  }
1150  }
1151  updated_bufs.emplace_back(&line_buf);
1152 }
1153 
1154 void VisualizationSceneVector3d::PrepareDisplacedMesh()
1155 {
1156  int dim = mesh->Dimension();
1157  int i, j, ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
1158  DenseMatrix pointmat;
1159  Array<int> vertices;
1160 
1161  // prepare the displaced mesh
1162  displine_buf.clear();
1163  gl3::GlBuilder builder = displine_buf.createBuilder();
1164 
1165  for (i = 0; i < ne; i++)
1166  {
1167  builder.glBegin(GL_LINE_LOOP);
1168  if (dim == 3)
1169  {
1170  mesh->GetBdrPointMatrix (i, pointmat);
1171  mesh->GetBdrElementVertices (i, vertices);
1172  }
1173  else
1174  {
1175  mesh->GetPointMatrix(i, pointmat);
1176  mesh->GetElementVertices(i, vertices);
1177  }
1178 
1179  for (j = 0; j < pointmat.Size(); j++)
1180  {
1181  pointmat(0, j) += (*solx)(vertices[j])*(ianimd)/ianimmax;
1182  pointmat(1, j) += (*soly)(vertices[j])*(ianimd)/ianimmax;
1183  pointmat(2, j) += (*solz)(vertices[j])*(ianimd)/ianimmax;
1184  }
1185 
1186  for (j = 0; j < pointmat.Size(); j++)
1187  {
1188  builder.glVertex3d (pointmat(0, j), pointmat(1, j), pointmat(2, j) );
1189  }
1190  builder.glEnd();
1191  }
1192  updated_bufs.emplace_back(&displine_buf);
1193 }
1194 
1195 void ArrowsDrawOrNot (Array<int> l[], int nv, Vector & sol,
1196  int nl, Array<double> & level)
1197 {
1198  static int first_time = 1;
1199  static int nll = nl;
1200 
1201  if (!first_time && nll == nl)
1202  {
1203  return;
1204  }
1205  else
1206  {
1207  first_time = 1;
1208  nll = nl;
1209  }
1210 
1211  int i,j;
1212  double v;
1213  double eps = (level[1] - level[0])/2;
1214 
1215  for (i = 0; i <= nl; i++)
1216  {
1217  l[i].SetSize(0);
1218  }
1219 
1220  for (j = 0; j < nv; j++)
1221  {
1222  v = sol(j);
1223  for (i = 0; i <= nl; i++)
1224  {
1225  if (fabs(v-level[i]) < eps)
1226  {
1227  l[i].Append(j);
1228  break;
1229  }
1230  if (v < level[i] - eps)
1231  {
1232  break;
1233  }
1234  }
1235  }
1236 }
1237 
1238 int ArrowDrawOrNot (double v, int nl, Array<double> & level)
1239 {
1240  double eps = (level[nl] - level[0])/10;
1241  for (int i = 0; i <= nl; i++)
1242  {
1243  if (fabs(v-level[i]) < eps)
1244  {
1245  return 1;
1246  }
1247  if (v < level[i] - eps)
1248  {
1249  return 0;
1250  }
1251  }
1252  return 0;
1253 }
1254 
1255 void VisualizationSceneVector3d::DrawVector(gl3::GlBuilder& builder,
1256  int type, double v0, double v1,
1257  double v2, double sx, double sy,
1258  double sz, double s)
1259 {
1260  static int nv = mesh -> GetNV();
1261  static double volume = (bb.x[1]-bb.x[0])*(bb.y[1]-bb.y[0])*(bb.z[1]-bb.z[0]);
1262  static double h = pow(volume/nv, 0.333);
1263  static double hh = pow(volume, 0.333) / 10;
1264 
1265  switch (type)
1266  {
1267  case 1:
1268  {
1269  arrow_type = 0;
1270  arrow_scaling_type = 0;
1271  // glColor3f(0, 0, 0); // color is set in Draw()
1272  Arrow(builder,v0,v1,v2,sx,sy,sz,s);
1273  }
1274  break;
1275 
1276  case 2:
1277  {
1278  arrow_type = 1;
1279  arrow_scaling_type = 1;
1280  MySetColor(builder, s, minv, maxv);
1281  Arrow(builder,v0,v1,v2,sx,sy,sz,h,0.125);
1282  }
1283  break;
1284 
1285  case 3:
1286  {
1287  arrow_type = 1;
1288  arrow_scaling_type = 1;
1289  // MySetColor(s,maxv,minv);
1290  MySetColor(builder, s, minv, maxv);
1291  Arrow(builder,v0,v1,v2,sx,sy,sz,h*s/maxv,0.125);
1292  }
1293  break;
1294 
1295  case 4:
1296  case 5:
1297  {
1298  arrow_type = 1;
1299  arrow_scaling_type = 1;
1300  builder.glColor3f(0.3, 0.3, 0.3);
1301  Arrow(builder,v0,v1,v2,sx,sy,sz,hh*s/maxv,0.125);
1302  }
1303  break;
1304  }
1305 }
1306 
1307 void VisualizationSceneVector3d::PrepareVectorField()
1308 {
1309  int i, nv = mesh -> GetNV();
1310  double *vertex;
1311 
1312  vector_buf.clear();
1313  gl3::GlBuilder builder = vector_buf.createBuilder();
1314 
1315  switch (drawvector)
1316  {
1317  case 0:
1318  break;
1319 
1320  case 1:
1321  for (i = 0; i < nv; i++)
1322  if (drawmesh != 2 || ArrowDrawOrNot((*sol)(i), nl, level))
1323  {
1324  vertex = mesh->GetVertex(i);
1325  DrawVector(builder, drawvector, vertex[0], vertex[1], vertex[2],
1326  (*solx)(i), (*soly)(i), (*solz)(i), (*sol)(i));
1327  }
1328  break;
1329 
1330  case 2:
1331  {
1332  arrow_type = 1;
1333  arrow_scaling_type = 1;
1334  for (i = 0; i < nv; i++)
1335  if (drawmesh != 2 || ArrowDrawOrNot((*sol)(i), nl, level))
1336  {
1337  vertex = mesh->GetVertex(i);
1338  DrawVector(builder, drawvector, vertex[0], vertex[1], vertex[2],
1339  (*solx)(i), (*soly)(i), (*solz)(i), (*sol)(i));
1340  }
1341  }
1342  break;
1343 
1344  case 3:
1345  {
1346  arrow_type = 1;
1347  arrow_scaling_type = 1;
1348 
1349  for (i = 0; i < nv; i++)
1350  if (drawmesh != 2 || ArrowDrawOrNot((*sol)(i), nl, level))
1351  {
1352  vertex = mesh->GetVertex(i);
1353  DrawVector(builder, drawvector, vertex[0], vertex[1], vertex[2],
1354  (*solx)(i), (*soly)(i), (*solz)(i), (*sol)(i));
1355  }
1356  }
1357  break;
1358 
1359  case 4:
1360  {
1361  Array<int> *l = new Array<int>[nl+1];
1362  ArrowsDrawOrNot(l, nv, *sol, nl, level);
1363 
1364  int j,k;
1365 
1366  for (k = 0; k < vflevel.Size(); k++)
1367  {
1368  i = vflevel[k];
1369  for (j = 0; j < l[i].Size(); j++)
1370  {
1371  vertex = mesh->GetVertex( l[i][j] );
1372  DrawVector(builder, drawvector, vertex[0], vertex[1], vertex[2],
1373  (*solx)(l[i][j]), (*soly)(l[i][j]), (*solz)(l[i][j]),
1374  (*sol)(l[i][j]));
1375  }
1376  }
1377 
1378  delete [] l;
1379  }
1380  break;
1381 
1382  case 5:
1383  {
1384  int dim = mesh->Dimension();
1385  int j, k, ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE();
1386  Array<int> vertices;
1387  Array<bool> vert_marker(nv);
1388 
1389  vert_marker = false;
1390  for (k = 0; k < ne; k++)
1391  {
1392  if (dim == 3)
1393  {
1394  mesh->GetBdrElementVertices(k, vertices);
1395  }
1396  else
1397  {
1398  mesh->GetElementVertices(k, vertices);
1399  }
1400  for (j = 0; j < vertices.Size(); j++)
1401  {
1402  i = vertices[j];
1403  if (vert_marker[i]) { continue; }
1404  vertex = mesh->GetVertex(i);
1405  DrawVector(builder, drawvector, vertex[0], vertex[1], vertex[2],
1406  (*solx)(i), (*soly)(i), (*solz)(i), (*sol)(i));
1407  vert_marker[i] = true;
1408  }
1409  }
1410  }
1411  break;
1412  }
1413  updated_bufs.emplace_back(&vector_buf);
1414 }
1415 
1416 void VisualizationSceneVector3d::PrepareCuttingPlane()
1417 {
1418  if (cp_drawelems == 0 || cplane != 1 || drawvector == 0 ||
1419  mesh->Dimension() != 3)
1420  {
1421  VisualizationSceneSolution3d::PrepareCuttingPlane();
1422  return;
1423  }
1424 
1425  int i, j, k, n = 0;
1426  int flag[4], ind[6][2]= {{0,3},{0,2},{0,1},{1,2},{1,3},{2,3}};
1427  double t, point[4][4], val[4][3];
1428 
1429  cplane_buf.clear();
1430  gl3::GlBuilder builder = cplane_buf.createBuilder();
1431 
1432  DenseMatrix pointmat(3,4);
1433  double * coord;
1434 
1435  Array<int> nodes;
1436  for (i = 0; i < mesh->GetNE(); i++)
1437  {
1438  if (mesh->GetElementType(i) != Element::TETRAHEDRON)
1439  {
1440  continue;
1441  }
1442  // TODO: support for hex elements as in
1443  // VisualizationSceneSolution3d::CuttingPlaneFunc
1444 
1445  n = 0; // n will be the number of intersection points
1446  mesh -> GetElementVertices(i,nodes);
1447  mesh -> GetPointMatrix (i,pointmat);
1448  for (j=0; j<4; j++)
1449  if (node_pos[nodes[j]] == 0.0)
1450  {
1451  flag[j] = 1;
1452  coord = mesh -> GetVertex(nodes[j]);
1453  for (k=0; k<3; k++)
1454  {
1455  point[n][k] = coord[k];
1456  }
1457 
1458  point[n][3] = (*sol)(nodes[j]);
1459  val[n][0] = (*solx)(nodes[j]);
1460  val[n][1] = (*soly)(nodes[j]);
1461  val[n][2] = (*solz)(nodes[j]);
1462  n++;
1463  }
1464  else if (node_pos[nodes[j]] < 0.)
1465  {
1466  flag[j] = -1;
1467  }
1468  else
1469  {
1470  flag[j] = 0;
1471  }
1472 
1473  for (j=0; j<6; j++)
1474  if (flag[ind[j][0]] != 1 && flag[ind[j][1]] != 1 &&
1475  flag[ind[j][0]] != flag[ind[j][1]])
1476  {
1477  t = node_pos[ nodes[ind[j][1]] ] / (node_pos[ nodes[ind[j][1]] ] -
1478  node_pos[ nodes[ind[j][0]] ] );
1479  for (k=0; k<3; k++)
1480  point[n][k] = t*(pointmat(k,ind[j][0]) -
1481  pointmat(k,ind[j][1])) +
1482  pointmat(k,ind[j][1]);
1483 
1484  point[n][3] = t * ((*sol)(nodes[ind[j][0]]) -
1485  (*sol)(nodes[ind[j][1]])) +
1486  (*sol)(nodes[ind[j][1]]);
1487  val[n][0] = t * ((*solx)(nodes[ind[j][0]]) -
1488  (*solx)(nodes[ind[j][1]])) +
1489  (*solx)(nodes[ind[j][1]]);
1490  val[n][1] = t * ((*soly)(nodes[ind[j][0]]) -
1491  (*soly)(nodes[ind[j][1]])) +
1492  (*soly)(nodes[ind[j][1]]);
1493  val[n][2] = t * ((*solz)(nodes[ind[j][0]]) -
1494  (*solz)(nodes[ind[j][1]])) +
1495  (*solz)(nodes[ind[j][1]]);
1496  n++;
1497  }
1498 
1499  if (n > 2)
1500  {
1501  double v10[] = { point[1][0]-point[0][0],
1502  point[1][1]-point[0][1],
1503  point[1][2]-point[0][2]
1504  };
1505  double v21[] = { point[2][0]-point[1][0],
1506  point[2][1]-point[1][1],
1507  point[2][2]-point[1][2]
1508  };
1509 
1510  double norm[] = { v10[1]*v21[2]-v10[2]*v21[1],
1511  v10[2]*v21[0]-v10[0]*v21[2],
1512  v10[0]*v21[1]-v10[1]*v21[0]
1513  };
1514 
1515  double * eq = CuttingPlane -> Equation();
1516 
1517  if ( eq[0]*norm[0]+eq[1]*norm[1]+eq[2]*norm[2] > 0.0 )
1518  {
1519  if (drawvector != 5)
1520  {
1521  builder.glBegin (GL_POLYGON);
1522  for (j=0; j<n; j++)
1523  {
1524  MySetColor ( builder, point[j][3], maxv, minv);
1525  builder.glNormal3dv (CuttingPlane -> Equation());
1526  builder.glVertex3d (point[j][0],point[j][1],point[j][2]);
1527  }
1528  }
1529  builder.glEnd();
1530  if (drawvector)
1531  for (j=0; j<n; j++)
1532  DrawVector(builder, drawvector, point[n][0], point[n][1], point[n][2],
1533  val[n][0],val[n][1], val[n][2], point[n][3]);
1534  }
1535  else
1536  {
1537  if (drawvector != 5)
1538  {
1539  builder.glBegin (GL_POLYGON);
1540  for (j=n-1; j>=0; j--)
1541  {
1542  MySetColor ( builder, point[j][3], minv, maxv);
1543  builder.glNormal3dv (CuttingPlane -> Equation());
1544  builder.glVertex3d (point[j][0],point[j][1],point[j][2]);
1545  }
1546  }
1547  builder.glEnd();
1548  if (drawvector)
1549  for (j=n-1; j>=0; j--)
1550  DrawVector(builder, drawvector, point[n][0], point[n][1], point[n][2],
1551  val[n][0],val[n][1], val[n][2], point[n][3]);
1552  }
1553  }
1554  }
1555  updated_bufs.emplace_back(&cplane_buf);
1556 }
1557 
1558 gl3::SceneInfo VisualizationSceneVector3d::GetSceneObjs()
1559 {
1560  if (colorbar)
1561  {
1562  Array<double>* cb_level = nullptr;
1563  if (drawvector == 4)
1564  {
1565  cb_level = &dvflevel;
1566  }
1567  else if (drawmesh == 2 || cp_drawmesh >= 2)
1568  {
1569  cb_level = &level;
1570  }
1571  PrepareColorBar(minv, maxv, cb_level);
1572  }
1573  gl3::SceneInfo scene = VisualizationSceneScalarData::GetSceneObjs();
1574  gl3::RenderParams params = GetMeshDrawParams();
1575  params.use_clip_plane = cplane;
1576  double* cp_eqn = CuttingPlane->Equation();
1577  params.clip_plane_eqn = {cp_eqn[0], cp_eqn[1], cp_eqn[2], cp_eqn[3]};
1578  params.contains_translucent = matAlpha < 1.0;
1579  // draw vector field
1580  if (drawvector == 2 || drawvector == 3)
1581  {
1582  scene.queue.emplace_back(params, &vector_buf);
1583  }
1584  // draw elements
1585  if (drawelems)
1586  {
1587  scene.queue.emplace_back(params, &disp_buf);
1588  }
1589 
1590  if (cplane && cp_drawelems)
1591  {
1592  params.use_clip_plane = false;
1593  scene.queue.emplace_back(params, &cplane_buf);
1594  params.use_clip_plane = true;
1595  }
1596 
1597  params.contains_translucent = false;
1598  params.num_pt_lights = 0;
1599 
1600  if (drawvector > 3)
1601  {
1602  scene.queue.emplace_back(params, &vector_buf);
1603  }
1604 
1605  params.mesh_material = VisualizationScene::BLK_MAT;
1606  params.static_color = GetLineColor();
1607 
1608  if (drawvector == 1)
1609  {
1610  scene.queue.emplace_back(params, &vector_buf);
1611  }
1612 
1613  // draw lines
1614  if (drawmesh)
1615  {
1616  scene.queue.emplace_back(params, &line_buf);
1617  }
1618 
1619  // draw displacement
1620  if (drawdisp)
1621  {
1622  params.static_color = {1.f, 0.f, 0.f, 1.f};
1623  scene.queue.emplace_back(params, &displine_buf);
1624  params.static_color = GetLineColor();
1625  }
1626  // ruler may have mixture of polygons and lines
1627  if (cp_drawmesh)
1628  {
1629  params.use_clip_plane = false;
1630  scene.queue.emplace_back(params, &cplines_buf);
1631  }
1632  return scene;
1633 }
SendExposeEvent
void SendExposeEvent()
Send expose event. In our case MyReshape is executed and Draw after it.
Definition: aux_vis.cpp:346
KeyvPressed
void KeyvPressed()
Definition: vsvector.cpp:141
GLVisGeometryRefiner
thread_local GeometryRefiner GLVisGeometryRefiner
Definition: glvis.cpp:71
VisualizationSceneVector3d::DrawVector
void DrawVector(gl3::GlBuilder &builder, int type, double v0, double v1, double v2, double sx, double sy, double sz, double s)
Definition: vsvector3d.cpp:1255
VisualizationSceneVector3d::PrepareCuttingPlane
virtual void PrepareCuttingPlane()
Definition: vsvector3d.cpp:1416
VisualizationSceneVector3d::VisualizationSceneVector3d
VisualizationSceneVector3d(Mesh &m, Vector &sx, Vector &sy, Vector &sz)
Definition: vsvector3d.cpp:341
gl3::GlBuilder::glVertex3d
void glVertex3d(double x, double y, double z)
Definition: types.hpp:332
Normalize
int Normalize(double v[])
Definition: geom_utils.hpp:38
cut_TriGeoms
thread_local Array< int > cut_TriGeoms
Definition: vssolution3d.cpp:34
KeyrPressed
void KeyrPressed()
Definition: vsdata.cpp:531
wnd
thread_local SdlWindow * wnd
Definition: aux_vis.cpp:50
gl3::GlBuilder::glColor3f
void glColor3f(float r, float g, float b)
Definition: types.hpp:425
cut_QuadPts
thread_local IntegrationRule cut_QuadPts
Definition: vssolution3d.cpp:31
gl3::GlBuilder::glEnd
void glEnd()
Definition: types.hpp:309
VisualizationSceneScalarData::GetSceneObjs
virtual gl3::SceneInfo GetSceneObjs()
Definition: vsdata.cpp:984
VisualizationSceneVector3d::PrepareFlat
virtual void PrepareFlat()
Definition: vsvector3d.cpp:518
VisualizationSceneVector3d::PrepareLines
virtual void PrepareLines()
Definition: vsvector3d.cpp:898
ianimmax
thread_local int ianimmax
Definition: vsvector.cpp:108
ianim
thread_local int ianim
Definition: vsvector.cpp:107
VisualizationSceneVector3d::GetHelpString
virtual std::string GetHelpString() const
Definition: vsvector3d.cpp:33
gl3::GlBuilder::glNormal3dv
void glNormal3dv(const double *d)
Definition: types.hpp:407
KeyNPressed
void KeyNPressed()
Definition: vsvector.cpp:116
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
CutReferenceElements
void CutReferenceElements(int TimesToRefine, double lambda)
Definition: vssolution3d.cpp:1609
KeyDPressed
void KeyDPressed()
Definition: vsvector.cpp:110
VisualizationSceneVector3d::GetSceneObjs
virtual gl3::SceneInfo GetSceneObjs()
Definition: vsvector3d.cpp:1558
ArrowDrawOrNot
int ArrowDrawOrNot(double v, int nl, Array< double > &level)
Definition: vsvector3d.cpp:1238
gl3::GlBuilder::glNormal3d
void glNormal3d(double nx, double ny, double nz)
Definition: types.hpp:401
KeyRPressed
void KeyRPressed()
Definition: vsdata.cpp:546
cut_QuadGeoms
thread_local Array< int > cut_QuadGeoms
Definition: vssolution3d.cpp:32
RemoveIdleFunc
void RemoveIdleFunc(void(*Func)(void))
Definition: aux_vis.cpp:497
gl3::GlBuilder
Crude fixed-function OpenGL emulation helper.
Definition: types.hpp:261
extra_caption
thread_local string extra_caption
Definition: glvis.cpp:61
MainLoop
void MainLoop()
Definition: aux_vis.cpp:515
VisualizationSceneVector3d::NewMeshAndSolution
void NewMeshAndSolution(Mesh *new_m, GridFunction *new_v)
Definition: vsvector3d.cpp:453
vsvector3d
thread_local VisualizationSceneVector3d * vsvector3d
Definition: vsvector3d.cpp:114
KeyUPressed
void KeyUPressed()
Definition: vsvector.cpp:203
cut_TriPts
thread_local IntegrationRule cut_TriPts
Definition: vssolution3d.cpp:33
ArrowsDrawOrNot
void ArrowsDrawOrNot(Array< int > l[], int nv, Vector &sol, int nl, Array< double > &level)
Definition: vsvector3d.cpp:1195
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
VisualizationSceneVector3d::~VisualizationSceneVector3d
virtual ~VisualizationSceneVector3d()
Definition: vsvector3d.cpp:439
VisualizationSceneVector3d::PrepareVectorField
virtual void PrepareVectorField()
Definition: vsvector3d.cpp:1307
KeyuPressed
void KeyuPressed()
Definition: vsvector.cpp:158
gl3::RenderParams::static_color
std::array< float, 4 > static_color
Definition: renderer.hpp:48
SdlWindow::setOnKeyDown
void setOnKeyDown(int key, Delegate func)
Definition: sdl.hpp:188
gl3::SceneInfo::queue
RenderQueue queue
Definition: renderer.hpp:63
gl3::RenderParams::mesh_material
Material mesh_material
Definition: renderer.hpp:44
VisualizationSceneVector3d::ToggleVectorFieldLevel
void ToggleVectorFieldLevel(int v)
Definition: vsvector3d.cpp:204
gl3::RenderParams::clip_plane_eqn
std::array< double, 4 > clip_plane_eqn
Definition: renderer.hpp:52
locscene
thread_local VisualizationScene * locscene
Definition: aux_vis.cpp:38
KeyBPressed
void KeyBPressed()
Definition: vsvector.cpp:122
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