1 #include <mystdlib.h>
2
3 #include <myadt.hpp>
4 #include <meshing.hpp>
5 #include <csg.hpp>
6 #include <stlgeom.hpp>
7
8
9 // #include <parallel.hpp>
10 #include <visual.hpp>
11 #include <limits>
12 namespace netgen
13 {
14
15
GetVSSolution()16 DLL_HEADER VisualSceneSolution & GetVSSolution()
17 {
18 static VisualSceneSolution vssolution;
19 return vssolution;
20 }
21
22
23 // extern shared_ptr<Mesh> mesh;
24 extern VisualSceneMesh vsmesh;
25
26
AddUserVisualizationObject(UserVisualizationObject * vis)27 DLL_HEADER void AddUserVisualizationObject (UserVisualizationObject * vis)
28 {
29 // vssolution.AddUserVisualizationObject (vis);
30 GetVSSolution().AddUserVisualizationObject (vis);
31 }
DeleteUserVisualizationObject(UserVisualizationObject * vis)32 DLL_HEADER void DeleteUserVisualizationObject (UserVisualizationObject * vis)
33 {
34 // vssolution.AddUserVisualizationObject (vis);
35 GetVSSolution().DeleteUserVisualizationObject (vis);
36 }
37
38
SolData()39 VisualSceneSolution :: SolData :: SolData ()
40 : data (0), solclass(0)
41 { ; }
42
~SolData()43 VisualSceneSolution :: SolData :: ~SolData ()
44 {
45 // delete [] name;
46 delete data;
47 delete solclass;
48 }
49
50
VisualSceneSolution()51 VisualSceneSolution :: VisualSceneSolution ()
52 : VisualScene()
53 {
54 // cout << "init VisualSceneSolution" << endl;
55 surfellist = 0;
56 linelist = 0;
57 element1dlist = 0;
58 clipplanelist_scal = 0;
59 clipplanelist_vec = 0;
60 isolinelist = 0;
61 clipplane_isolinelist = 0;
62 surface_vector_list = 0;
63 isosurface_list = 0;
64
65 fieldlineslist = 0;
66 pointcurvelist = 0;
67
68 num_fieldlineslists = 0;
69
70
71 surfeltimestamp = GetTimeStamp();
72 surfellinetimestamp = GetTimeStamp();
73 clipplanetimestamp = GetTimeStamp();
74 solutiontimestamp = GetTimeStamp();
75 fieldlinestimestamp = GetTimeStamp();
76 pointcurve_timestamp = GetTimeStamp();
77 surface_vector_timestamp = GetTimeStamp();
78 isosurface_timestamp = GetTimeStamp();
79 timetimestamp = GetTimeStamp();
80 // AddVisualizationScene ("solution", &vssolution);
81 }
82
~VisualSceneSolution()83 VisualSceneSolution :: ~VisualSceneSolution ()
84 {
85 // cout << "exit VisualSceneSolution" << endl;
86 ClearSolutionData();
87 }
88
89 /*
90 void VisualSceneSolution :: SetMesh (shared_ptr<Mesh> amesh)
91 {
92 wp_mesh = amesh;
93 }
94 */
95
AddSolutionData(SolData * sd)96 void VisualSceneSolution :: AddSolutionData (SolData * sd)
97 {
98 shared_ptr<Mesh> mesh = GetMesh();
99
100 NgLock meshlock1 (mesh->MajorMutex(), 1);
101 int funcnr = -1;
102 for (int i = 0; i < soldata.Size(); i++)
103 {
104 // if (strcmp (soldata[i]->name, sd->name) == 0)
105 if (soldata[i]->name == sd->name)
106 {
107 delete soldata[i];
108 soldata[i] = sd;
109 funcnr = i;
110 break;
111 }
112 }
113
114 if (funcnr == -1)
115 {
116 soldata.Append (sd);
117 funcnr = soldata.Size()-1;
118 }
119
120 SolData * nsd = soldata[funcnr];
121
122 nsd->size = 0;
123 if (mesh)
124 {
125 switch (nsd->soltype)
126 {
127 case SOL_NODAL: nsd->size = mesh->GetNV(); break;
128 case SOL_ELEMENT: nsd->size = mesh->GetNE(); break;
129 case SOL_SURFACE_ELEMENT: nsd->size = mesh->GetNSE(); break;
130 case SOL_NONCONTINUOUS:
131 {
132 switch (nsd->order)
133 {
134 case 0: nsd->size = mesh->GetNE(); break;
135 case 1: nsd->size = 6 * mesh->GetNE(); break;
136 case 2: nsd->size = 18 * mesh->GetNE(); break;
137 }
138 break;
139 }
140 case SOL_SURFACE_NONCONTINUOUS:
141 {
142 switch (nsd->order)
143 {
144 case 0: nsd->size = mesh->GetNSE(); break;
145 case 1: nsd->size = 4 * mesh->GetNSE(); break;
146 case 2: nsd->size = 9 * mesh->GetNSE(); break;
147 }
148 break;
149 }
150 default:
151 nsd->size = 0;
152 }
153 solutiontimestamp = NextTimeStamp();
154 }
155 }
156
157
ClearSolutionData()158 void VisualSceneSolution :: ClearSolutionData ()
159 {
160 for (int i = 0; i < soldata.Size(); i++)
161 delete soldata[i];
162 soldata.SetSize (0);
163 }
164
UpdateSolutionTimeStamp()165 void VisualSceneSolution :: UpdateSolutionTimeStamp ()
166 {
167 solutiontimestamp = NextTimeStamp();
168 }
169
GetSolData(int i)170 VisualSceneSolution::SolData * VisualSceneSolution :: GetSolData (int i)
171 {
172 if (i >= 0 && i < soldata.Size())
173 return soldata[i];
174 else
175 return NULL;
176 }
177
178
179
180
SaveSolutionData(const char * filename)181 void VisualSceneSolution :: SaveSolutionData (const char * filename)
182 {
183 shared_ptr<Mesh> mesh = GetMesh();
184
185 PrintMessage (1, "Write solution data to file ", filename);
186
187
188 if (strcmp (&filename[strlen(filename)-3], "sol") == 0)
189 {
190 ofstream ost(filename);
191 for (int i = 0; i < soldata.Size(); i++)
192 {
193 const SolData & sol = *soldata[i];
194
195 ost << "solution "
196 << sol.name
197 << " -size=" << sol.size
198 << " -components=" << sol.components
199 << " -order=" << sol.order;
200 if (sol.iscomplex)
201 ost << " -complex";
202
203 switch (sol.soltype)
204 {
205 case SOL_NODAL:
206 ost << " -type=nodal"; break;
207 case SOL_ELEMENT:
208 ost << " -type=element"; break;
209 case SOL_SURFACE_ELEMENT:
210 ost << " -type=surfaceelement"; break;
211 case SOL_NONCONTINUOUS:
212 ost << " -type=noncontinuous"; break;
213 case SOL_SURFACE_NONCONTINUOUS:
214 ost << " -type=surfacenoncontinuous"; break;
215 default:
216 cerr << "save solution data, case not handeld" << endl;
217 }
218
219 ost << endl;
220 for (int j = 0; j < sol.size; j++)
221 {
222 for (int k = 0; k < sol.components; k++)
223 ost << sol.data[j*sol.dist+k] << " ";
224 ost << "\n";
225 }
226 }
227 }
228
229
230 if (strcmp (&filename[strlen(filename)-3], "vtk") == 0)
231 {
232 string surf_fn = filename;
233 surf_fn.erase (strlen(filename)-4);
234 surf_fn += "_surf.vtk";
235
236 cout << "surface mesh = " << surf_fn << endl;
237
238 ofstream surf_ost(surf_fn.c_str());
239
240 surf_ost << "# vtk DataFile Version 1.0\n"
241 << "NGSolve surface mesh\n"
242 << "ASCII\n"
243 << "DATASET UNSTRUCTURED_GRID\n\n";
244
245 surf_ost << "POINTS " << mesh->GetNP() << " float\n";
246 for (PointIndex pi = PointIndex::BASE; pi < mesh->GetNP()+PointIndex::BASE; pi++)
247 {
248 const MeshPoint & mp = (*mesh)[pi];
249 surf_ost << mp(0) << " " << mp(1) << " " << mp(2) << "\n";
250 }
251
252 int cntverts = 0;
253 for (SurfaceElementIndex sei = 0; sei < mesh->GetNSE(); sei++)
254 cntverts += 1 + (*mesh)[sei].GetNP();
255
256 surf_ost << "\nCELLS " << mesh->GetNSE() << " " << cntverts << "\n";
257 for (SurfaceElementIndex sei = 0; sei < mesh->GetNSE(); sei++)
258 {
259 const Element2d & el = (*mesh)[sei];
260 surf_ost << el.GetNP();
261 for (int j = 0; j < el.GetNP(); j++)
262 surf_ost << " " << el[j] - PointIndex::BASE;
263 surf_ost << "\n";
264 }
265 surf_ost << "\nCELL_TYPES " << mesh->GetNSE() << "\n";
266 for (SurfaceElementIndex sei = 0; sei < mesh->GetNSE(); sei++)
267 {
268 const Element2d & el = (*mesh)[sei];
269 switch (el.GetType())
270 {
271 case QUAD: surf_ost << 9; break;
272 case TRIG: surf_ost << 5; break;
273 default:
274 cerr << "not implemented 2378" << endl;
275 }
276 surf_ost << "\n";
277 }
278
279
280
281 ofstream ost(filename);
282
283 ost << "# vtk DataFile Version 1.0\n"
284 << "NGSolve solution\n"
285 << "ASCII\n"
286 << "DATASET UNSTRUCTURED_GRID\n\n";
287
288 ost << "POINTS " << mesh->GetNP() << " float\n";
289 for (PointIndex pi = PointIndex::BASE; pi < mesh->GetNP()+PointIndex::BASE; pi++)
290 {
291 const MeshPoint & mp = (*mesh)[pi];
292 ost << mp(0) << " " << mp(1) << " " << mp(2) << "\n";
293 }
294
295 cntverts = 0;
296 for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
297 cntverts += 1 + (*mesh)[ei].GetNP();
298
299 ost << "\nCELLS " << mesh->GetNE() << " " << cntverts << "\n";
300 for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
301 {
302 const Element & el = (*mesh)[ei];
303 ost << el.GetNP();
304 for (int j = 0; j < el.GetNP(); j++)
305 ost << " " << el[j] - PointIndex::BASE;
306 ost << "\n";
307 }
308 ost << "\nCELL_TYPES " << mesh->GetNE() << "\n";
309 for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
310 {
311 const Element & el = (*mesh)[ei];
312 switch (el.GetType())
313 {
314 case TET: ost << 10; break;
315 default:
316 cerr << "not implemented 67324" << endl;
317 }
318 ost << "\n";
319 }
320
321
322 ost << "CELL_DATA " << mesh->GetNE() << "\n";
323 for (int i = 0; i < soldata.Size(); i++)
324 {
325 ost << "VECTORS bfield float\n";
326 SolutionData & sol = *(soldata[i] -> solclass);
327 double values[3];
328
329 for (int elnr = 0; elnr < mesh->GetNE(); elnr++)
330 {
331 sol.GetValue (elnr, 0.25, 0.25, 0.25, values);
332 ost << values[0] << " " << values[1] << " " << values[2] << "\n";
333 }
334 }
335
336 /*
337 ost << "POINT_DATA " << mesh->GetNP() << "\n";
338 for (int i = 0; i < soldata.Size(); i++)
339 {
340 ost << "VECTORS bfield float\n";
341 SolutionData & sol = *(soldata[i] -> solclass);
342
343 for (PointIndex pi = PointIndex::BASE;
344 pi < mesh->GetNP()+PointIndex::BASE; pi++)
345 {
346 double values[3], sumvalues[3] = { 0, 0, 0 };
347
348 NgFlatArray<int> els = mesh->GetTopology().GetVertexElements(pi);
349
350 for (int j = 0; j < els.Size(); j++)
351 {
352 sol.GetValue (els[j]-1, 0.25, 0.25, 0.25, values);
353 for (int k = 0; k < 3; k++)
354 sumvalues[k] += values[k];
355 }
356 for (int k = 0; k < 3; k++)
357 sumvalues[k] /= els.Size();
358
359 ost << sumvalues[0] << " " << sumvalues[1] << " " << sumvalues[2] << "\n";
360 }
361 }
362 */
363 }
364
365 }
366
367
368
369
DrawScene()370 void VisualSceneSolution :: DrawScene ()
371 {
372 try
373 {
374 shared_ptr<Mesh> mesh = GetMesh();
375
376 if (!mesh)
377 {
378 VisualScene::DrawScene();
379 return;
380 }
381
382 // static NgLock mem_lock(mem_mutex);
383 // mem_lock.Lock();
384
385 NgLock meshlock1 (mesh->MajorMutex(), true);
386 NgLock meshlock (mesh->Mutex(), true);
387
388 BuildScene();
389
390 CreateTexture (numtexturecols, lineartexture, 0.5, GL_MODULATE);
391
392 glClearColor(backcolor, backcolor, backcolor, 1);
393 // glClearColor(backcolor, backcolor, backcolor, 0);
394 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
395
396 SetLight();
397
398 glPushMatrix();
399 glMultMatrixd (transformationmat);
400
401 glMatrixMode (GL_MODELVIEW);
402
403 glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
404
405 glPolygonOffset (1, 1);
406
407 glEnable (GL_POLYGON_OFFSET_FILL);
408
409 glEnable (GL_COLOR_MATERIAL);
410
411 if (usetexture)
412 {
413 SetTextureMode (usetexture);
414
415 glMatrixMode (GL_TEXTURE);
416 glLoadIdentity();
417
418 if (usetexture == 1)
419 {
420 double hmax = maxval;
421 double hmin = minval;
422 if (invcolor) Swap (hmax, hmin);
423
424 if (fabs (hmax - hmin) > 1e-30)
425 glScaled (1.0 / (hmin - hmax), 0, 0);
426 else
427 glScaled (1e30, 0, 0);
428
429 glTranslatef (-hmax, 0, 0);
430 }
431 else
432 {
433 glTranslatef (0.5, 0, 0);
434 glRotatef(360 * netgen::GetVSSolution().time, 0, 0, -1);
435 if (fabs (maxval) > 1e-10)
436 glScalef(0.5/maxval, 0.5/maxval, 0.5/maxval);
437 else
438 glScalef (1e10, 1e10, 1e10);
439 }
440 glMatrixMode (GL_MODELVIEW);
441 }
442
443
444
445
446 if (vispar.drawfilledtrigs || vispar.drawtetsdomain > 0 || vispar.drawdomainsurf > 0)
447 {
448 // Change for Martin:
449
450 // orig:
451 SetClippingPlane ();
452
453 glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
454 // glEnable(GL_BLEND);
455 glDisable(GL_BLEND);
456 glCallList (surfellist);
457
458 #ifdef USE_BUFFERS
459 // static int timer = NgProfiler::CreateTimer ("Solution::drawing - DrawSurfaceElements VBO");
460 // NgProfiler::StartTimer(timer);
461 glEnableClientState(GL_VERTEX_ARRAY);
462 glEnableClientState(GL_NORMAL_ARRAY);
463 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
464 glDrawElements(GL_TRIANGLES, surfel_vbo_size, GL_UNSIGNED_INT, 0);
465 glDisableClientState(GL_VERTEX_ARRAY);
466 glDisableClientState(GL_NORMAL_ARRAY);
467 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
468 // NgProfiler::StopTimer(timer);
469 #endif
470
471 glDisable(GL_BLEND);
472 /*
473 // transparent test ...
474 glColor4f (1, 0, 0, 0.1);
475 glEnable (GL_COLOR_MATERIAL);
476
477 glDepthFunc(GL_GREATER);
478 glDepthMask(GL_FALSE);
479 // glBlendFunc(GL_ONE_MINUS_DST_ALPHA,GL_DST_ALPHA);
480 glBlendFunc(GL_ONE_MINUS_SRC_ALPHA,GL_SRC_ALPHA);
481
482 glCallList (surfellist);
483
484 glDisable(GL_BLEND);
485 glDepthFunc(GL_LEQUAL);
486 glDepthMask(GL_TRUE);
487
488 glCallList (surfellist);
489 // end test ...
490 */
491
492
493 glCallList (surface_vector_list);
494 glDisable(GL_CLIP_PLANE0);
495 }
496
497
498 if (showclipsolution)
499 {
500 if (clipsolution == 1)
501 {
502 // Martin
503 // orig:
504 glCallList (clipplanelist_scal);
505
506 // transparent experiments
507 // see http://wiki.delphigl.com/index.php/Blenden
508
509 /*
510 glColor4f (1, 1, 1, 0.5);
511 glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
512 glEnable(GL_BLEND);
513 glEnable(GL_COLOR);
514 glDepthFunc(GL_GREATER);
515 glDepthMask(GL_FALSE);
516
517 glCallList (clipplanelist_scal);
518 glDepthFunc(GL_LEQUAL);
519 glDepthMask(GL_TRUE);
520
521 glCallList (clipplanelist_scal);
522 glDisable(GL_BLEND);
523 */
524
525
526 /*
527 // latest transparent version ...
528 glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
529 glEnable(GL_BLEND);
530 glEnable(GL_DEPTH_TEST);
531
532 // CreateTexture (numtexturecols, lineartexture, 0.25, GL_MODULATE);
533 // glCallList (clipplanelist_scal);
534
535 glEnable(GL_BLEND);
536 // glDisable(GL_DEPTH_TEST);
537
538 // CreateTexture (numtexturecols, lineartexture, 0.25, GL_MODULATE);
539 glCallList (clipplanelist_scal);
540
541
542 // glDepthFunc(GL_LEQUAL);
543 // glDepthMask(GL_TRUE);
544 // glCallList (clipplanelist_scal);
545 glEnable(GL_DEPTH_TEST);
546 glDisable(GL_BLEND);
547 */
548 // end test
549 }
550 if (clipsolution == 2)
551 {
552 // glDisable(GL_DEPTH_TEST);
553 glCallList (clipplanelist_vec);
554 // glEnable(GL_DEPTH_TEST);
555 }
556 }
557
558
559
560 if (draw_fieldlines)
561 {
562 SetClippingPlane();
563 if (num_fieldlineslists <= 1)
564 glCallList (fieldlineslist);
565 else
566 { // animated
567 int start = int (time / 10 * num_fieldlineslists);
568 for (int ln = 0; ln < 10; ln++)
569 {
570 int nr = fieldlineslist + (start + ln) % num_fieldlineslists;
571 glCallList (nr);
572 }
573 }
574 glDisable(GL_CLIP_PLANE0);
575 }
576
577 if(drawpointcurves)
578 {
579 glCallList(pointcurvelist);
580 }
581
582
583 glMatrixMode (GL_TEXTURE);
584 glLoadIdentity();
585 glMatrixMode (GL_MODELVIEW);
586
587 glDisable (GL_TEXTURE_1D);
588 glDisable (GL_TEXTURE_2D);
589
590 glDisable (GL_POLYGON_OFFSET_FILL);
591 glDisable (GL_COLOR_MATERIAL);
592
593 if (draw_isosurface)
594 glCallList (isosurface_list);
595
596
597 GLfloat matcol0[] = { 0, 0, 0, 1 };
598 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcol0);
599 glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, matcol0);
600 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, matcol0);
601
602 glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
603 glLineWidth (1.0f);
604 glColor3f (0.0f, 0.0f, 0.0f);
605 glDisable (GL_LINE_SMOOTH);
606
607 if (vispar.drawedges)
608 {
609 glCallList (element1dlist);
610 }
611
612 if (vispar.drawoutline && !numisolines)
613 {
614 SetClippingPlane ();
615 glDepthMask(GL_FALSE);
616 glCallList (linelist);
617 glDepthMask(GL_TRUE);
618
619 glDisable(GL_CLIP_PLANE0);
620 }
621
622 if (numisolines)
623 {
624 SetClippingPlane ();
625 glCallList (isolinelist);
626
627 glDisable(GL_CLIP_PLANE0);
628 glCallList (clipplane_isolinelist);
629 }
630
631
632 // user visualization
633
634 for (int i = 0; i < user_vis.Size(); i++)
635 user_vis[i] -> Draw();
636
637 glPopMatrix();
638
639 glDisable(GL_CLIP_PLANE0);
640 DrawColorBar (minval, maxval, logscale, lineartexture);
641
642 if (vispar.drawcoordinatecross)
643 DrawCoordinateCross ();
644 DrawNetgenLogo ();
645
646 glFinish();
647
648
649 // delete lock;
650 // mem_lock.UnLock();
651 }
652 catch (bad_weak_ptr e)
653 {
654 // cout << "don't have a mesh to visualize" << endl;
655 VisualScene::DrawScene();
656 }
657 }
658
659
660
661 /*
662 void VisualSceneSolution :: RealVec3d (const double * values, Vec3d & v,
663 bool iscomplex, bool imag)
664 {
665 if (!iscomplex)
666 {
667 v.X() = values[0];
668 v.Y() = values[1];
669 v.Z() = values[2];
670 }
671 else
672 {
673 if (!imag)
674 {
675 v.X() = values[0];
676 v.Y() = values[2];
677 v.Z() = values[4];
678 }
679 else
680 {
681 v.X() = values[1];
682 v.Y() = values[3];
683 v.Z() = values[5];
684 }
685 }
686 }
687 */
RealVec3d(const double * values,bool iscomplex,bool imag)688 Vec<3> VisualSceneSolution :: RealVec3d (const double * values,
689 bool iscomplex, bool imag)
690 {
691 Vec<3> v;
692 if (!iscomplex)
693 {
694 for (int j = 0; j < 3; j++)
695 v(j) = values[j];
696 }
697 else
698 {
699 if (!imag)
700 {
701 for (int j = 0; j < 3; j++)
702 v(j) = values[2*j];
703 }
704 else
705 {
706 for (int j = 0; j < 3; j++)
707 v(j) = values[2*j+1];
708 }
709 }
710 return v;
711 }
712
713
RealVec3d(const double * values,Vec3d & v,bool iscomplex,double phaser,double phasei)714 void VisualSceneSolution :: RealVec3d (const double * values, Vec3d & v,
715 bool iscomplex, double phaser, double phasei)
716 {
717 if (!iscomplex)
718 {
719 v.X() = values[0];
720 v.Y() = values[1];
721 v.Z() = values[2];
722 }
723 else
724 {
725 for (int i = 0; i < 3; i++)
726 v.X(i+1) = phaser * values[2*i] + phasei * values[2*i+1];
727 }
728 }
729
730
731
732
BuildScene(int zoomall)733 void VisualSceneSolution :: BuildScene (int zoomall)
734 {
735 try
736 {
737 shared_ptr<Mesh> mesh = GetMesh();
738
739 if (!mesh)
740 {
741 VisualScene::BuildScene (zoomall);
742 return;
743 }
744
745 /*
746 if (!cone_list)
747 {
748 cone_list = glGenLists (1);
749 glNewList (cone_list, GL_COMPILE);
750 DrawCone (Point<3> (0,0,0), Point<3> (0,0,1), 0.4);
751 glEndList();
752 }
753 */
754
755 // vispar.colormeshsize = 1;
756
757 // recalc clipping plane
758 SetClippingPlane ();
759 glDisable(GL_CLIP_PLANE0);
760
761
762 SolData * sol = NULL;
763 SolData * vsol = NULL;
764
765 if (scalfunction != -1)
766 sol = soldata[scalfunction];
767 if (vecfunction != -1)
768 vsol = soldata[vecfunction];
769
770 if (mesh->GetTimeStamp () > solutiontimestamp)
771 {
772 sol = NULL;
773 vsol = NULL;
774 }
775
776
777 if (sol && sol->solclass) sol->solclass->SetMultiDimComponent (multidimcomponent);
778 if (vsol && vsol->solclass) vsol->solclass->SetMultiDimComponent (multidimcomponent);
779
780 if (!autoscale || (!sol && !vsol) )
781 {
782 minval = mminval;
783 maxval = mmaxval;
784 }
785 else
786 {
787 if (mesh->GetTimeStamp () > surfeltimestamp ||
788 vispar.clipping.timestamp > clipplanetimestamp ||
789 solutiontimestamp > surfeltimestamp)
790 {
791 GetMinMax (scalfunction, scalcomp, minval, maxval);
792 }
793 }
794
795 if (mesh->GetTimeStamp() > surfeltimestamp ||
796 solutiontimestamp > surfeltimestamp ||
797 zoomall)
798 {
799 if (mesh->GetTimeStamp() > surfeltimestamp || zoomall)
800 {
801 // mesh has changed
802
803 Point3d pmin, pmax;
804 static double oldrad = 0;
805
806 mesh->GetBox (pmin, pmax, -1);
807 center = Center (pmin, pmax);
808 rad = 0.5 * Dist (pmin, pmax);
809
810 glEnable (GL_NORMALIZE);
811
812 if (rad > 1.5 * oldrad ||
813 mesh->GetMajorTimeStamp() > surfeltimestamp ||
814 zoomall)
815 {
816 CalcTransformationMatrices();
817 oldrad = rad;
818 }
819 }
820
821 DrawSurfaceElements();
822
823 surfeltimestamp = max2 (solutiontimestamp, mesh->GetTimeStamp());
824 }
825
826 if (mesh->GetTimeStamp() > surfellinetimestamp ||
827 subdivision_timestamp > surfellinetimestamp ||
828 (deform && solutiontimestamp > surfellinetimestamp) ||
829 zoomall)
830 {
831 DrawSurfaceElementLines();
832 surfellinetimestamp = max2 (solutiontimestamp, mesh->GetTimeStamp());
833 }
834
835
836 if (vispar.drawedges)
837 Draw1DElements();
838
839
840
841 if (mesh->GetTimeStamp() > surface_vector_timestamp ||
842 solutiontimestamp > surface_vector_timestamp ||
843 zoomall)
844 {
845 if (surface_vector_list)
846 glDeleteLists (surface_vector_list, 1);
847
848 surface_vector_list = glGenLists (1);
849 glNewList (surface_vector_list, GL_COMPILE);
850
851 glEnable (GL_NORMALIZE);
852 DrawSurfaceVectors();
853
854 glEndList ();
855
856 surface_vector_timestamp =
857 max2 (mesh->GetTimeStamp(), solutiontimestamp);
858 }
859
860 if (clipplanetimestamp < vispar.clipping.timestamp ||
861 clipplanetimestamp < solutiontimestamp)
862 {
863
864 // cout << "clipsolution = " << clipsolution << endl;
865 if (vispar.clipping.enable && clipsolution == 2)
866 {
867 mesh->Mutex().unlock();
868 mesh->BuildElementSearchTree();
869 mesh->Mutex().lock();
870 }
871
872
873 if (vispar.clipping.enable && clipsolution == 1 && sol)
874 DrawClipPlaneTrigs ();
875
876 if (clipplanelist_vec)
877 glDeleteLists (clipplanelist_vec, 1);
878
879 clipplanelist_vec = glGenLists (1);
880 glNewList (clipplanelist_vec, GL_COMPILE);
881
882 if (vispar.clipping.enable && clipsolution == 2 && vsol)
883 {
884 SetTextureMode (usetexture);
885
886 if (autoscale)
887 GetMinMax (vecfunction, 0, minval, maxval);
888
889 NgArray<ClipPlanePoint> cpp;
890 GetClippingPlaneGrid (cpp);
891
892 for (int i = 0; i < cpp.Size(); i++)
893 {
894 const ClipPlanePoint & p = cpp[i];
895 double values[6];
896 Vec3d v;
897
898 bool drawelem =
899 GetValues (vsol, p.elnr, p.lami(0), p.lami(1), p.lami(2), values);
900 // RealVec3d (values, v, vsol->iscomplex, imag_part);
901 v = RealVec3d (values, vsol->iscomplex, imag_part);
902
903 double val = v.Length();
904
905 if (drawelem && val > 1e-10 * maxval)
906 {
907 v *= (rad / val / gridsize * 0.5);
908
909 SetOpenGlColor (val);
910 DrawCone (p.p, p.p+v, rad / gridsize * 0.2);
911 }
912 }
913 }
914
915 glEndList ();
916 }
917
918
919 if (mesh->GetTimeStamp() > isosurface_timestamp ||
920 solutiontimestamp > isosurface_timestamp ||
921 zoomall)
922 {
923 if (isosurface_list)
924 glDeleteLists (isosurface_list, 1);
925
926 isosurface_list = glGenLists (1);
927 glNewList (isosurface_list, GL_COMPILE);
928
929 glEnable (GL_NORMALIZE);
930 DrawIsoSurface(sol, vsol, scalcomp);
931
932 glEndList ();
933
934 isosurface_timestamp =
935 max2 (mesh->GetTimeStamp(), solutiontimestamp);
936 }
937
938 if(mesh->GetTimeStamp() > pointcurve_timestamp ||
939 solutiontimestamp > pointcurve_timestamp)
940 {
941 if(pointcurvelist)
942 glDeleteLists(pointcurvelist,1);
943
944
945 if(mesh->GetNumPointCurves() > 0)
946 {
947 pointcurvelist = glGenLists(1);
948 glNewList(pointcurvelist,GL_COMPILE);
949 //glColor3f (1.0f, 0.f, 0.f);
950
951 for(int i=0; i<mesh->GetNumPointCurves(); i++)
952 {
953 Box3d box;
954 box.SetPoint(mesh->GetPointCurvePoint(i,0));
955 for(int j=1; j<mesh->GetNumPointsOfPointCurve(i); j++)
956 box.AddPoint(mesh->GetPointCurvePoint(i,j));
957 double diam = box.CalcDiam();
958
959 double thick = min2(0.1*diam, 0.001*rad);
960
961 double red,green,blue;
962 mesh->GetPointCurveColor(i,red,green,blue);
963 glColor3f (red, green, blue);
964 for(int j=0; j<mesh->GetNumPointsOfPointCurve(i)-1; j++)
965 {
966 DrawCylinder(mesh->GetPointCurvePoint(i,j),
967 mesh->GetPointCurvePoint(i,j+1),
968 thick);
969 }
970 }
971 glEndList();
972 }
973
974 }
975
976
977 if (
978 numisolines &&
979 (clipplanetimestamp < vispar.clipping.timestamp ||
980 clipplanetimestamp < solutiontimestamp)
981 )
982 {
983 if (isolinelist) glDeleteLists (isolinelist, 1);
984
985 isolinelist = glGenLists (1);
986 glNewList (isolinelist, GL_COMPILE);
987
988 Point<3> points[1100];
989 double values[1100];
990
991 int nse = mesh->GetNSE();
992
993 CurvedElements & curv = mesh->GetCurvedElements();
994
995 if (sol)
996 {
997 glBegin (GL_LINES);
998
999 for (SurfaceElementIndex sei = 0; sei < nse; sei++)
1000 {
1001 const Element2d & el = (*mesh)[sei];
1002
1003 bool curved = curv.IsHighOrder(); // && curv.IsSurfaceElementCurved(sei);
1004
1005 if (el.GetType() == TRIG || el.GetType() == TRIG6)
1006 {
1007 Point<3> lp1, lp2, lp3;
1008 if (!curved)
1009 {
1010 GetPointDeformation (el[0]-1, lp1);
1011 GetPointDeformation (el[1]-1, lp2);
1012 GetPointDeformation (el[2]-1, lp3);
1013 }
1014
1015 int n = 1 << subdivisions;
1016 int ii = 0;
1017 int ix, iy;
1018 for (iy = 0; iy <= n; iy++)
1019 for (ix = 0; ix <= n-iy; ix++)
1020 {
1021 double x = double(ix) / n;
1022 double y = double(iy) / n;
1023
1024 // TODO: consider return value (bool: draw/don't draw element)
1025 GetSurfValue (sol, sei, -1, x, y, scalcomp, values[ii]);
1026 Point<2> xref(x,y);
1027
1028 if (curved)
1029 mesh->GetCurvedElements().
1030 CalcSurfaceTransformation (xref, sei, points[ii]);
1031 else
1032 points[ii] = lp3 + x * (lp1-lp3) + y * (lp2-lp3);
1033
1034 if (deform)
1035 {
1036 points[ii] += GetSurfDeformation (sei, -1, x, y);
1037 }
1038 ii++;
1039 }
1040
1041 ii = 0;
1042 for (iy = 0; iy < n; iy++, ii++)
1043 for (ix = 0; ix < n-iy; ix++, ii++)
1044 {
1045 int index[] = { ii, ii+1, ii+n-iy+1,
1046 ii+1, ii+n-iy+2, ii+n-iy+1 };
1047
1048 DrawIsoLines (points[index[0]], points[index[1]], points[index[2]],
1049 values[index[0]], values[index[1]], values[index[2]]);
1050
1051 if (ix < n-iy-1)
1052 DrawIsoLines (points[index[3]], points[index[4]], points[index[5]],
1053 values[index[3]], values[index[4]], values[index[5]]);
1054 }
1055 }
1056
1057
1058 if (el.GetType() == QUAD || el.GetType() == QUAD6 || el.GetType() == QUAD8 )
1059 {
1060 Point<3> lpi[4];
1061 Vec<3> vx = 0.0, vy = 0.0, vtwist = 0.0, def;
1062 if (!curved)
1063 {
1064 for (int j = 0; j < 4; j++)
1065 GetPointDeformation (el[j]-1, lpi[j]);
1066 vx = lpi[1]-lpi[0];
1067 vy = lpi[3]-lpi[0];
1068 vtwist = (lpi[0]-lpi[1]) + (lpi[2]-lpi[3]);
1069 }
1070
1071 int n = 1 << subdivisions;
1072 int ix, iy, ii = 0;
1073 for (iy = 0; iy <= n; iy++)
1074 for (ix = 0; ix <= n; ix++, ii++)
1075 {
1076 double x = double(ix) / n;
1077 double y = double(iy) / n;
1078
1079 // TODO: consider return value (bool: draw/don't draw element)
1080 GetSurfValue (sol, sei, -1, x, y, scalcomp, values[ii]);
1081 Point<2> xref(x,y);
1082
1083 if (curved)
1084 mesh->GetCurvedElements().
1085 CalcSurfaceTransformation (xref, sei, points[ii]);
1086 else
1087 points[ii] = lpi[0] + x * vx + y * vy + x*y * vtwist;
1088
1089 if (deform)
1090 points[ii] += GetSurfDeformation (sei, -1, x, y);
1091 }
1092
1093 ii = 0;
1094 for (iy = 0; iy < n; iy++, ii++)
1095 for (ix = 0; ix < n; ix++, ii++)
1096 {
1097 DrawIsoLines (points[ii], points[ii+1], points[ii+n+1],
1098 values[ii], values[ii+1], values[ii+n+1]);
1099 DrawIsoLines (points[ii+1], points[ii+n+2], points[ii+n+1],
1100 values[ii+1], values[ii+n+2], values[ii+n+1]);
1101 }
1102 }
1103 }
1104 glEnd();
1105 }
1106 glEndList ();
1107
1108 if (clipplane_isolinelist) glDeleteLists (clipplane_isolinelist, 1);
1109
1110 if (vispar.clipping.enable && clipsolution == 1 && sol)
1111 {
1112 clipplane_isolinelist = glGenLists (1);
1113 glNewList (clipplane_isolinelist, GL_COMPILE);
1114
1115 NgArray<ClipPlaneTrig> cpt;
1116 NgArray<ClipPlanePoint> pts;
1117 GetClippingPlaneTrigs (cpt, pts);
1118 bool drawelem;
1119
1120 glNormal3d (-clipplane[0], -clipplane[1], -clipplane[2]);
1121 glBegin (GL_LINES);
1122
1123 if (numisolines)
1124 for (int i = 0; i < cpt.Size(); i++)
1125 {
1126 const ClipPlaneTrig & trig = cpt[i];
1127 double vali[3];
1128 for (int j = 0; j < 3; j++)
1129 {
1130 Point<3> lami = pts[trig.points[j].pnr].lami;
1131 drawelem = GetValue (sol, trig.elnr, lami(0), lami(1), lami(2),
1132 scalcomp, vali[j]);
1133 }
1134 if ( drawelem )
1135 DrawIsoLines (pts[trig.points[0].pnr].p,
1136 pts[trig.points[1].pnr].p,
1137 pts[trig.points[2].pnr].p,
1138 vali[0], vali[1], vali[2]);
1139 }
1140 glEnd();
1141 glEndList ();
1142 }
1143 }
1144
1145 clipplanetimestamp = max2 (vispar.clipping.timestamp, solutiontimestamp);
1146 }
1147 catch (bad_weak_ptr e)
1148 {
1149 PrintMessage (3, "vssolution::buildscene: don't have a mesh to visualize");
1150 VisualScene::BuildScene (zoomall);
1151 }
1152 }
1153
Draw1DElements()1154 void VisualSceneSolution :: Draw1DElements ()
1155 {
1156 shared_ptr<Mesh> mesh = GetMesh();
1157
1158 if (element1dlist)
1159 glDeleteLists (element1dlist, 1);
1160
1161 element1dlist = glGenLists (1);
1162 glNewList (element1dlist, GL_COMPILE);
1163
1164 int npt = (1 << subdivisions) + 1;
1165 NgArray<double> pref(npt), values(npt);
1166 NgArray<Point<3> > points(npt);
1167
1168 const SolData * sol = NULL;
1169 if (scalfunction != -1) sol = soldata[scalfunction];
1170
1171 const SolData * vsol = NULL;
1172 if (deform && vecfunction != -1) vsol = soldata[vecfunction];
1173
1174 int ncomp = 0;
1175 if (sol) ncomp = sol->components;
1176 if (vsol) ncomp = vsol->components;
1177 NgArray<double> mvalues(ncomp);
1178
1179
1180 for (int i = 0; i < npt; i++)
1181 pref[i] = double(i) / (npt-1);
1182 int meshdim = mesh->GetDimension();
1183 for (SegmentIndex i = 0; i < mesh -> GetNSeg(); i++)
1184 {
1185 // mesh->GetCurvedElements().
1186 // CalcMultiPointSegmentTransformation (&pref, i, &points, NULL);
1187 // const Segment & seg = mesh -> LineSegment(i);
1188 for (int j = 0; j < npt; j++)
1189 mesh->GetCurvedElements().
1190 CalcSegmentTransformation (pref[j], i, points[j]);
1191 if (vsol)
1192 {
1193 for (int j = 0; j < npt; j++)
1194 {
1195 vsol->solclass->GetSegmentValue (i, pref[j], &mvalues[0]);
1196 // values[j] = ExtractValue (sol, scalcomp, &mvalues[0]);
1197 for (int k = 0; k < min(ncomp, 3); k++)
1198 points[j](k) += scaledeform * mvalues[k];
1199 // points[j](0) += scaledeform * mvalues[0];
1200 // points[j](1) += scaledeform * mvalues[1];
1201 }
1202 }
1203 else if (sol)
1204 {
1205 for (int j = 0; j < npt; j++)
1206 {
1207 sol->solclass->GetSegmentValue (i, pref[j], &mvalues[0]);
1208 values[j] = ExtractValue (sol, scalcomp, &mvalues[0]);
1209 if (meshdim <= 2)
1210 points[j](meshdim) += scaledeform * values[j];
1211 }
1212 }
1213
1214 glBegin (GL_LINE_STRIP);
1215 for (int i = 0; i < npt; i++)
1216 glVertex3dv (points[i]);
1217 glEnd();
1218 }
1219
1220 glEndList ();
1221 }
1222
DrawSurfaceElements()1223 void VisualSceneSolution :: DrawSurfaceElements ()
1224 {
1225 shared_ptr<Mesh> mesh = GetMesh();
1226
1227 static int timer = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements");
1228 /*
1229 static int timerstart = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements start");
1230 static int timerloops = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements loops");
1231 static int timerlist = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements list");
1232 static int timerbuffer = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements buffer");
1233 static int timer1 = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1");
1234 static int timer1a = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1a");
1235 static int timer1b = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1b");
1236 static int timer1c = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1c");
1237 static int timer2 = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 2");
1238 static int timer2a = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 2a");
1239 static int timer2b = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 2b");
1240 */
1241 NgProfiler::RegionTimer reg (timer);
1242
1243
1244 #ifdef PARALLELGL
1245
1246 if (id == 0 && ntasks > 1)
1247 {
1248 InitParallelGL();
1249
1250 par_surfellists.SetSize (ntasks);
1251
1252 MyMPI_SendCmd ("redraw");
1253 MyMPI_SendCmd ("solsurfellist");
1254
1255 for ( int dest = 1; dest < ntasks; dest++ )
1256 MyMPI_Recv (par_surfellists[dest], dest, MPI_TAG_VIS);
1257
1258 if (surfellist)
1259 glDeleteLists (surfellist, 1);
1260
1261 surfellist = glGenLists (1);
1262 glNewList (surfellist, GL_COMPILE);
1263
1264 for ( int dest = 1; dest < ntasks; dest++ )
1265 glCallList (par_surfellists[dest]);
1266
1267 glEndList();
1268 return;
1269 }
1270 #endif
1271
1272 // NgProfiler::StartTimer(timerstart);
1273
1274 if (surfellist)
1275 glDeleteLists (surfellist, 1);
1276
1277 surfellist = glGenLists (1);
1278 glNewList (surfellist, GL_COMPILE);
1279
1280
1281 const SolData * sol = NULL;
1282
1283 if (scalfunction != -1)
1284 sol = soldata[scalfunction];
1285
1286 if (mesh->GetTimeStamp () > solutiontimestamp)
1287 sol = NULL;
1288
1289 if (sol && sol->solclass) sol->solclass->SetMultiDimComponent (multidimcomponent);
1290
1291
1292
1293 glLineWidth (1.0f);
1294
1295 GLfloat col_grey[] = { 0.6f, 0.6f, 0.6f, 1.0f };
1296 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, col_grey);
1297
1298
1299 int nse = mesh->GetNSE();
1300
1301 SetTextureMode (usetexture);
1302
1303 CurvedElements & curv = mesh->GetCurvedElements();
1304
1305 int n = 1 << subdivisions;
1306 int npt = sqr(n+1);
1307
1308 NgArray<Point<2> > pref (npt);
1309 NgArray<Point<3> > points (npt);
1310 NgArray<Mat<3,2> > dxdxis (npt);
1311 NgArray<Vec<3> > nvs(npt);
1312 NgArray<double> values(npt);
1313
1314 NgArray<double> mvalues(npt);
1315 int sol_comp = (sol && sol->draw_surface) ? sol->components : 0;
1316 NgArray<Point<2,SIMD<double>> > simd_pref ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
1317 NgArray<Point<3,SIMD<double>> > simd_points ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
1318 NgArray<Mat<3,2,SIMD<double>> > simd_dxdxis ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
1319 NgArray<Vec<3,SIMD<double>> > simd_nvs( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
1320 NgArray<SIMD<double>> simd_values( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() * sol_comp);
1321
1322
1323
1324 // NgArray<Point<3,float>> glob_pnts;
1325 // NgArray<Vec<3,float>> glob_nvs;
1326 // NgArray<double> glob_values;
1327
1328 if (sol && sol->draw_surface) mvalues.SetSize (npt * sol->components);
1329
1330 NgArray<complex<double> > valuesc(npt);
1331
1332 #ifdef USE_BUFFERS
1333 if (has_surfel_vbo)
1334 glDeleteBuffers (4, &surfel_vbo[0]);
1335 glGenBuffers (4, &surfel_vbo[0]);
1336
1337 has_surfel_vbo = true;
1338 glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[0]);
1339 glBufferData (GL_ARRAY_BUFFER,
1340 nse*npt*sizeof(Point<3,double>),
1341 NULL, GL_STATIC_DRAW);
1342 glVertexPointer(3, GL_DOUBLE, 0, 0);
1343 // glEnableClientState(GL_VERTEX_ARRAY);
1344
1345 glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[1]);
1346 glBufferData (GL_ARRAY_BUFFER,
1347 nse*npt*sizeof(Vec<3,double>),
1348 NULL, GL_STATIC_DRAW);
1349 // glEnableClientState(GL_NORMAL_ARRAY);
1350 glNormalPointer(GL_DOUBLE, 0, 0);
1351
1352 // glEnableClientState(GL_TEXTURE_COORD_ARRAY);
1353 glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[2]);
1354 glBufferData (GL_ARRAY_BUFFER, nse*npt*sizeof(double), NULL, GL_STATIC_DRAW);
1355 glTexCoordPointer(1, GL_DOUBLE, 0, 0);
1356
1357 glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]);
1358 glBufferData(GL_ELEMENT_ARRAY_BUFFER, nse*npt*6*sizeof(int), NULL, GL_STATIC_DRAW);
1359 surfel_vbo_size = 0;
1360 #endif
1361
1362
1363 // NgProfiler::StopTimer(timerstart);
1364
1365 for (SurfaceElementIndex sei = 0; sei < nse; sei++)
1366 {
1367 const Element2d & el = (*mesh)[sei];
1368
1369 if (vispar.drawdomainsurf > 0)
1370 {
1371 if (mesh->GetDimension() == 3)
1372 {
1373 if (vispar.drawdomainsurf != mesh->GetFaceDescriptor(el.GetIndex()).DomainIn() &&
1374 vispar.drawdomainsurf != mesh->GetFaceDescriptor(el.GetIndex()).DomainOut())
1375 continue;
1376 }
1377 else
1378 {
1379 if (el.GetIndex() != vispar.drawdomainsurf) continue;
1380 }
1381 }
1382
1383
1384
1385 if ( el.GetType() == QUAD || el.GetType() == QUAD6 || el.GetType() == QUAD8 )
1386 {
1387 bool curved = curv.IsSurfaceElementCurved (sei);
1388
1389
1390 for (int iy = 0, ii = 0; iy <= n; iy++)
1391 for (int ix = 0; ix <= n; ix++, ii++)
1392 pref[ii] = Point<2> (double(ix)/n, double(iy)/n);
1393
1394 int npt = (n+1)*(n+1);
1395 if (curved)
1396 {
1397 for (int ii = 0; ii < npt; ii++)
1398 {
1399 Point<2> xref = pref[ii];
1400
1401 mesh->GetCurvedElements().
1402 CalcSurfaceTransformation (xref, sei, points[ii], dxdxis[ii]);
1403 nvs[ii] = Cross (dxdxis[ii].Col(0), dxdxis[ii].Col(1));
1404 nvs[ii].Normalize();
1405 }
1406 }
1407 else
1408 {
1409 Point<3> lpi[4];
1410 Vec<3> vx, vy, vtwist;
1411
1412 for (int k = 0; k < 4; k++)
1413 GetPointDeformation (el[k]-1, lpi[k]);
1414
1415 vx = lpi[1]-lpi[0];
1416 vy = lpi[3]-lpi[0];
1417 vtwist = (lpi[0]-lpi[1]) + (lpi[2]-lpi[3]);
1418
1419 for (int ii = 0; ii < npt; ii++)
1420 {
1421 double x = pref[ii](0);
1422 double y = pref[ii](1);
1423 points[ii] = lpi[0] + x * vx + y * vy + x*y * vtwist;
1424 for (int j = 0; j < 3; j++)
1425 {
1426 dxdxis[ii](j,0) = vx(j) + y*vtwist(j);
1427 dxdxis[ii](j,1) = vy(j) + x*vtwist(j);
1428 }
1429 }
1430
1431 Vec<3> nv = Cross (vx, vy);
1432 nv.Normalize();
1433 for (int ii = 0; ii < npt; ii++)
1434 nvs[ii] = nv;
1435 }
1436
1437
1438 bool drawelem = false;
1439 /*
1440 if (sol && sol->draw_surface)
1441 {
1442 if (usetexture == 2)
1443 for (int ii = 0; ii < npt; ii++)
1444 drawelem = GetSurfValueComplex (sol, sei, -1, pref[ii](0), pref[ii](1), scalcomp, valuesc[ii]);
1445 else
1446 for (int ii = 0; ii < npt; ii++)
1447 drawelem = GetSurfValue (sol, sei, -1, pref[ii](0), pref[ii](1), scalcomp, values[ii]);
1448 }
1449 */
1450 if (sol && sol->draw_surface)
1451 {
1452 drawelem = GetMultiSurfValues (sol, sei, -1, npt,
1453 &pref[0](0), &pref[1](0)-&pref[0](0),
1454 &points[0](0), &points[1](0)-&points[0](0),
1455 &dxdxis[0](0), &dxdxis[1](0)-&dxdxis[0](0),
1456 &mvalues[0], sol->components);
1457 if (usetexture == 2)
1458 for (int ii = 0; ii < npt; ii++)
1459 valuesc[ii] = ExtractValueComplex(sol, scalcomp, &mvalues[ii*sol->components]);
1460 else
1461 for (int ii = 0; ii < npt; ii++)
1462 values[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]);
1463 }
1464
1465
1466 if (deform)
1467 for (int ii = 0; ii < npt; ii++)
1468 points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
1469
1470
1471 int save_usetexture = usetexture;
1472 if (!drawelem)
1473 {
1474 usetexture = 0;
1475 SetTextureMode (0);
1476 }
1477
1478 int ii = 0;
1479
1480 glBegin (GL_QUADS);
1481
1482 for (int iy = 0; iy < n; iy++, ii++)
1483 for (int ix = 0; ix < n; ix++, ii++)
1484 {
1485 int index[] = { ii, ii+1, ii+n+2, ii+n+1 };
1486
1487 for (int j = 0; j < 4; j++)
1488 {
1489 if (drawelem)
1490 {
1491 if (usetexture != 2)
1492 SetOpenGlColor (values[index[j]]);
1493 else
1494 glTexCoord2f ( valuesc[index[j]].real(),
1495 valuesc[index[j]].imag() );
1496 }
1497 else
1498 glColor4fv (col_grey);
1499
1500 glNormal3dv (nvs[index[j]]);
1501 glVertex3dv (points[index[j]]);
1502 }
1503 }
1504 glEnd();
1505
1506 if (!drawelem && (usetexture != save_usetexture))
1507 {
1508 usetexture = save_usetexture;
1509 SetTextureMode (usetexture);
1510 }
1511
1512 }
1513 }
1514
1515 n = 1 << subdivisions;
1516 double invn = 1.0 / n;
1517 npt = (n+1)*(n+2)/2;
1518 // NgProfiler::StartTimer(timerloops);
1519 size_t base_pi = 0;
1520
1521 for (int iy = 0, ii = 0; iy <= n; iy++)
1522 for (int ix = 0; ix <= n-iy; ix++, ii++)
1523 pref[ii] = Point<2> (ix*invn, iy*invn);
1524
1525 constexpr size_t simd_size = SIMD<double>::Size();
1526 size_t simd_npt = (npt+simd_size-1)/simd_size;
1527
1528 for (size_t i = 0; i < simd_npt; i++)
1529 {
1530 simd_pref[i](0) = [&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](0) : 0; };
1531 simd_pref[i](1) = [&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](1) : 0; };
1532 }
1533
1534 NgArray<int> ind_reftrig;
1535 for (int iy = 0, ii = 0; iy < n; iy++,ii++)
1536 for (int ix = 0; ix < n-iy; ix++, ii++)
1537 {
1538 int nv = (ix+iy+1 < n) ? 6 : 3;
1539 int ind[] = { ii, ii+1, ii+n-iy+1,
1540 ii+n-iy+1, ii+1, ii+n-iy+2 };
1541 for (int j = 0; j < nv; j++)
1542 ind_reftrig.Append (ind[j]);
1543 }
1544 NgArray<int> glob_ind;
1545 glob_ind.SetSize(ind_reftrig.Size());
1546
1547
1548 for(SurfaceElementIndex sei = 0; sei < nse; sei++)
1549 {
1550 const Element2d & el = (*mesh)[sei];
1551 // if (el.GetIndex() <= 1) continue;
1552
1553 if(vispar.drawdomainsurf > 0)
1554 {
1555 if (mesh->GetDimension() == 3)
1556 {
1557 if (vispar.drawdomainsurf != mesh->GetFaceDescriptor(el.GetIndex()).DomainIn() &&
1558 vispar.drawdomainsurf != mesh->GetFaceDescriptor(el.GetIndex()).DomainOut())
1559 continue;
1560 }
1561 else
1562 {
1563 if (el.GetIndex() != vispar.drawdomainsurf)
1564 continue;
1565 }
1566 }
1567
1568 if ( el.GetType() == TRIG || el.GetType() == TRIG6 )
1569 {
1570 // NgProfiler::StartTimer(timer1);
1571 #ifdef __AVX_try_it_out__
1572 // NgProfiler::StartTimer(timer1a);
1573 bool curved = curv.IsSurfaceElementCurved(sei);
1574
1575 if (curved)
1576 {
1577 mesh->GetCurvedElements().
1578 CalcMultiPointSurfaceTransformation<3> (sei, simd_npt,
1579 &simd_pref[0](0), 2,
1580 &simd_points[0](0), 3,
1581 &simd_dxdxis[0](0,0), 6);
1582
1583 for (size_t ii = 0; ii < simd_npt; ii++)
1584 simd_nvs[ii] = Cross (simd_dxdxis[ii].Col(0), simd_dxdxis[ii].Col(1)).Normalize();
1585 }
1586 else
1587 {
1588 Point<3,SIMD<double>> p1 = mesh->Point (el[0]);
1589 Point<3,SIMD<double>> p2 = mesh->Point (el[1]);
1590 Point<3,SIMD<double>> p3 = mesh->Point (el[2]);
1591
1592 Vec<3,SIMD<double>> vx = p1-p3;
1593 Vec<3,SIMD<double>> vy = p2-p3;
1594 for (size_t ii = 0; ii < simd_npt; ii++)
1595 {
1596 simd_points[ii] = p3 + simd_pref[ii](0) * vx + simd_pref[ii](1) * vy;
1597 for (size_t j = 0; j < 3; j++)
1598 {
1599 simd_dxdxis[ii](j,0) = vx(j);
1600 simd_dxdxis[ii](j,1) = vy(j);
1601 }
1602 }
1603
1604 Vec<3,SIMD<double>> nv = Cross (vx, vy).Normalize();
1605 for (size_t ii = 0; ii < simd_npt; ii++)
1606 simd_nvs[ii] = nv;
1607 }
1608
1609
1610 bool drawelem = false;
1611 if (sol && sol->draw_surface)
1612 {
1613 // NgProfiler::StopTimer(timer1a);
1614 // NgProfiler::StartTimer(timer1b);
1615 drawelem = sol->solclass->GetMultiSurfValue (sei, -1, simd_npt,
1616 &simd_pref[0](0).Data(),
1617 &simd_points[0](0).Data(),
1618 &simd_dxdxis[0](0).Data(),
1619 &simd_values[0].Data());
1620 // NgProfiler::StopTimer(timer1b);
1621 // NgProfiler::StartTimer(timer1c);
1622
1623 for (size_t j = 0; j < sol->components; j++)
1624 for (size_t i = 0; i < npt; i++)
1625 mvalues[i*sol->components+j] = ((double*)&simd_values[j*simd_npt])[i];
1626
1627 if (usetexture == 2)
1628 for (int ii = 0; ii < npt; ii++)
1629 valuesc[ii] = ExtractValueComplex(sol, scalcomp, &mvalues[ii*sol->components]);
1630 else
1631 for (int ii = 0; ii < npt; ii++)
1632 values[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]);
1633 }
1634
1635 for (size_t i = 0; i < npt; i++)
1636 {
1637 size_t ii = i/4;
1638 size_t r = i%4;
1639 for (int j = 0; j < 2; j++)
1640 pref[i](j) = simd_pref[ii](j)[r];
1641 for (int j = 0; j < 3; j++)
1642 points[i](j) = simd_points[ii](j)[r];
1643 for (int j = 0; j < 3; j++)
1644 nvs[i](j) = simd_nvs[ii](j)[r];
1645 }
1646
1647 if (deform)
1648 for (int ii = 0; ii < npt; ii++)
1649 points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
1650
1651 // NgProfiler::StopTimer(timer1c);
1652
1653 #else
1654 bool curved = (*mesh)[sei].IsCurved();
1655
1656 for (int iy = 0, ii = 0; iy <= n; iy++)
1657 for (int ix = 0; ix <= n-iy; ix++, ii++)
1658 pref[ii] = Point<2> (ix*invn, iy*invn);
1659
1660 if (curved)
1661 {
1662 mesh->GetCurvedElements().
1663 CalcMultiPointSurfaceTransformation (&pref, sei, &points, &dxdxis);
1664
1665 for (int ii = 0; ii < npt; ii++)
1666 nvs[ii] = Cross (dxdxis[ii].Col(0), dxdxis[ii].Col(1)).Normalize();
1667 }
1668 else
1669 {
1670 Point<3> p1 = mesh->Point (el[0]);
1671 Point<3> p2 = mesh->Point (el[1]);
1672 Point<3> p3 = mesh->Point (el[2]);
1673
1674 Vec<3> vx = p1-p3;
1675 Vec<3> vy = p2-p3;
1676 for (int ii = 0; ii < npt; ii++)
1677 {
1678 points[ii] = p3 + pref[ii](0) * vx + pref[ii](1) * vy;
1679 for (int j = 0; j < 3; j++)
1680 {
1681 dxdxis[ii](j,0) = vx(j);
1682 dxdxis[ii](j,1) = vy(j);
1683 }
1684 }
1685
1686 Vec<3> nv = Cross (vx, vy).Normalize();
1687 for (int ii = 0; ii < npt; ii++)
1688 nvs[ii] = nv;
1689 }
1690
1691 bool drawelem = false;
1692 if (sol && sol->draw_surface)
1693 {
1694 drawelem = GetMultiSurfValues (sol, sei, -1, npt,
1695 &pref[0](0), &pref[1](0)-&pref[0](0),
1696 &points[0](0), &points[1](0)-&points[0](0),
1697 &dxdxis[0](0), &dxdxis[1](0)-&dxdxis[0](0),
1698 &mvalues[0], sol->components);
1699 if (usetexture == 2)
1700 for (int ii = 0; ii < npt; ii++)
1701 valuesc[ii] = ExtractValueComplex(sol, scalcomp, &mvalues[ii*sol->components]);
1702 else
1703 for (int ii = 0; ii < npt; ii++)
1704 values[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]);
1705 }
1706
1707 if (deform)
1708 for (int ii = 0; ii < npt; ii++)
1709 points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
1710 #endif
1711 // NgProfiler::StopTimer(timer1);
1712
1713 int save_usetexture = usetexture;
1714 if (!drawelem)
1715 {
1716 usetexture = 0;
1717 SetTextureMode (usetexture);
1718 }
1719
1720 // NgProfiler::StartTimer(timer2);
1721
1722 #ifdef USE_BUFFERS
1723 if (drawelem && usetexture == 1 && !logscale)
1724 {
1725 glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[0]);
1726 glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(Point<3,double>),
1727 npt*sizeof(Point<3,double>), &points[0][0]);
1728 glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[1]);
1729 glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(Vec<3,double>),
1730 npt*sizeof(Vec<3,double>), &nvs[0][0]);
1731 glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[2]);
1732 glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(double),
1733 npt*sizeof(double), &values[0]);
1734
1735 for (size_t i = 0; i < ind_reftrig.Size(); i++)
1736 glob_ind[i] = base_pi+ind_reftrig[i];
1737
1738 glBindBuffer (GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]);
1739 glBufferSubData (GL_ELEMENT_ARRAY_BUFFER, surfel_vbo_size*sizeof(int),
1740 ind_reftrig.Size()*sizeof(int), &glob_ind[0]);
1741 surfel_vbo_size += ind_reftrig.Size();
1742 base_pi += npt;
1743 }
1744
1745 else
1746 #endif
1747 for (int iy = 0, ii = 0; iy < n; iy++)
1748 {
1749 glBegin (GL_TRIANGLE_STRIP);
1750 for (int ix = 0; ix <= n-iy; ix++, ii++)
1751 for (int k = 0; k < 2; k++)
1752 {
1753 if (ix+iy+k > n) continue;
1754 int hi = (k == 0) ? ii : ii+n-iy+1;
1755 if (drawelem)
1756 {
1757 if (usetexture != 2)
1758 SetOpenGlColor (values[hi]);
1759 else
1760 glTexCoord2f ( valuesc[hi].real(), valuesc[hi].imag() );
1761 }
1762 else
1763 glColor4fv (col_grey);
1764
1765 glNormal3dv (nvs[hi]);
1766 glVertex3dv (points[hi]);
1767 }
1768 glEnd();
1769 }
1770
1771 // NgProfiler::StopTimer(timer2);
1772
1773
1774
1775
1776 if (!drawelem && (usetexture != save_usetexture))
1777 {
1778 usetexture = save_usetexture;
1779 SetTextureMode (usetexture);
1780 }
1781 }
1782 }
1783 // NgProfiler::StopTimer(timerloops);
1784
1785 // NgProfiler::StartTimer(timerbuffer);
1786
1787 // glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]);
1788 // glBufferData(GL_ELEMENT_ARRAY_BUFFER, glob_ind.Size()*sizeof(int), &glob_ind[0], GL_STATIC_DRAW);
1789 // surfel_vbo_size = glob_ind.Size();
1790
1791 // NgProfiler::StopTimer(timerbuffer);
1792
1793 // glDrawElements(GL_TRIANGLES, surfel_vbo_size, GL_UNSIGNED_INT, 0);
1794
1795 // glDrawElements(GL_TRIANGLES, glob_ind.Size(), GL_UNSIGNED_INT, &glob_ind[0]);
1796
1797 // glDisableClientState(GL_VERTEX_ARRAY);
1798 // glDisableClientState(GL_NORMAL_ARRAY);
1799 // glDisableClientState(GL_TEXTURE_COORD_ARRAY);
1800
1801 // glDeleteBuffers (1, &IndexVBOID);
1802 // glDeleteBuffers (4, &vboId[0]);
1803
1804
1805 // NgProfiler::StartTimer(timerlist);
1806 glEndList ();
1807 // NgProfiler::StopTimer(timerlist);
1808
1809 #ifdef PARALLELGL
1810 glFinish();
1811 if (id > 0)
1812 MyMPI_Send (surfellist, 0, MPI_TAG_VIS);
1813 #endif
1814 }
1815
1816
DrawSurfaceElementLines()1817 void VisualSceneSolution :: DrawSurfaceElementLines ()
1818 {
1819 shared_ptr<Mesh> mesh = GetMesh();
1820
1821 #ifdef PARALLELGL
1822 if (id == 0 && ntasks > 1)
1823 {
1824 InitParallelGL();
1825
1826 par_surfellists.SetSize (ntasks);
1827
1828 MyMPI_SendCmd ("redraw");
1829 MyMPI_SendCmd ("solsurfellinelist");
1830
1831 for ( int dest = 1; dest < ntasks; dest++ )
1832 MyMPI_Recv (par_surfellists[dest], dest, MPI_TAG_VIS);
1833
1834 if (linelist)
1835 glDeleteLists (linelist, 1);
1836
1837 linelist = glGenLists (1);
1838 glNewList (linelist, GL_COMPILE);
1839
1840 for ( int dest = 1; dest < ntasks; dest++ )
1841 glCallList (par_surfellists[dest]);
1842
1843 glEndList();
1844 return;
1845 }
1846 #endif
1847
1848 if (linelist)
1849 glDeleteLists (linelist, 1);
1850
1851 linelist = glGenLists (1);
1852 glNewList (linelist, GL_COMPILE);
1853
1854 glLineWidth (1.0f);
1855
1856 int nse = mesh->GetNSE();
1857 CurvedElements & curv = mesh->GetCurvedElements();
1858
1859 int n = 1 << subdivisions;
1860 NgArrayMem<Point<2>, 65> ptsloc(n+1);
1861 NgArrayMem<Point<3>, 65> ptsglob(n+1);
1862
1863 double trigpts[3][2] = { { 0, 0 }, { 0, 1 }, { 1, 0} };
1864 double trigvecs[3][2] = { { 1, 0 }, { 0, -1 }, { -1, 1} };
1865
1866 double quadpts[4][2] = { { 0, 0 }, { 1, 1 }, { 0, 1}, { 1, 0 } };
1867 double quadvecs[4][2] = { { 1, 0 }, { -1, 0}, { 0, -1}, { 0, 1 } };
1868
1869 for (SurfaceElementIndex sei = 0; sei < nse; sei++)
1870 {
1871 Element2d & el = (*mesh)[sei];
1872
1873 int nv = (el.GetType() == TRIG || el.GetType() == TRIG6) ? 3 : 4;
1874 for (int k = 0; k < nv; k++)
1875 {
1876 Point<2> p0;
1877 Vec<2> vtau;
1878 if (nv == 3)
1879 {
1880 p0 = Point<2>(trigpts[k][0], trigpts[k][1]);
1881 vtau = Vec<2>(trigvecs[k][0], trigvecs[k][1]);
1882 }
1883 else
1884 {
1885 p0 = Point<2>(quadpts[k][0], quadpts[k][1]);
1886 vtau = Vec<2>(quadvecs[k][0], quadvecs[k][1]);
1887 }
1888
1889 glBegin (GL_LINE_STRIP);
1890
1891 for (int ix = 0; ix <= n; ix++)
1892 ptsloc[ix] = p0 + (double(ix) / n) * vtau;
1893
1894 curv.CalcMultiPointSurfaceTransformation (&ptsloc, sei, &ptsglob, 0);
1895
1896 for (int ix = 0; ix <= n; ix++)
1897 {
1898 if (deform)
1899 ptsglob[ix] += GetSurfDeformation (sei, k, ptsloc[ix](0), ptsloc[ix](1));
1900 glVertex3dv (ptsglob[ix]);
1901 }
1902
1903 glEnd ();
1904 }
1905 }
1906 glEndList ();
1907
1908
1909 #ifdef PARALLELGL
1910 glFinish();
1911 if (id > 0)
1912 MyMPI_Send (linelist, 0, MPI_TAG_VIS);
1913 #endif
1914 }
1915
1916
1917
1918
1919
1920
1921
1922
1923
DrawIsoSurface(const SolData * sol,const SolData * vsol,int comp)1924 void VisualSceneSolution :: DrawIsoSurface(const SolData * sol,
1925 const SolData * vsol,
1926 int comp)
1927 {
1928 shared_ptr<Mesh> mesh = GetMesh();
1929
1930 if (!draw_isosurface) return;
1931 if (!sol) return;
1932
1933
1934 SetTextureMode (0);
1935 glColor3d (1.0, 0, 0);
1936 glEnable (GL_COLOR_MATERIAL);
1937
1938
1939 glBegin (GL_TRIANGLES);
1940
1941 int ne = mesh->GetNE();
1942
1943 const int edgei[6][2] =
1944 { { 0, 1 }, { 0, 2 }, { 0, 3 },
1945 { 1, 2 }, { 1, 3 }, { 2, 3 } };
1946
1947 double edgelam[6];
1948 Point<3> edgep[6];
1949 Vec<3> normp[6];
1950 double nodevali[4];
1951
1952 int cntce;
1953 int cpe1 = 0, cpe2 = 0, cpe3 = 0;
1954
1955 int n = 1 << subdivisions;
1956 int n3 = (n+1)*(n+1)*(n+1);
1957
1958 NgArray<Point<3> > grid(n3);
1959 NgArray<Point<3> > locgrid(n3);
1960 NgArray<Mat<3,3> > trans(n3);
1961 NgArray<double> val1(n3*sol->components);
1962 NgArray<Vec<3> > grads1(n3);
1963 NgArray<int> compress(n3);
1964
1965 MatrixFixWidth<3> pointmat(8);
1966 grads1 = Vec<3> (0.0);
1967
1968 for (ElementIndex ei = 0; ei < ne; ei++)
1969 {
1970 // if(vispar.clipdomain > 0 && vispar.clipdomain != (*mesh)[ei].GetIndex()) continue;
1971 // if(vispar.donotclipdomain > 0 && vispar.donotclipdomain == (*mesh)[ei].GetIndex()) continue;
1972
1973 ELEMENT_TYPE type = (*mesh)[ei].GetType();
1974 if (type == HEX || type == PRISM || type == TET || type == PYRAMID)
1975 {
1976 const Element & el = (*mesh)[ei];
1977
1978 int ii = 0;
1979 int cnt_valid = 0;
1980
1981 for (int ix = 0; ix <= n; ix++)
1982 for (int iy = 0; iy <= n; iy++)
1983 for (int iz = 0; iz <= n; iz++, ii++)
1984 {
1985 Point<3> ploc;
1986 compress[ii] = ii;
1987
1988 switch (type)
1989 {
1990 case PRISM:
1991 if (ix+iy <= n)
1992 {
1993 ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
1994 compress[ii] = cnt_valid;
1995 cnt_valid++;
1996 }
1997 else
1998 compress[ii] = -1;
1999 break;
2000 case TET:
2001 if (ix+iy+iz <= n)
2002 {
2003 ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
2004 compress[ii] = cnt_valid;
2005 cnt_valid++;
2006 }
2007 else
2008 compress[ii] = -1;
2009 break;
2010 case HEX:
2011 ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
2012 break;
2013 case PYRAMID:
2014 ploc = Point<3> (double(ix) / n * (1-double(iz)/n),
2015 double(iy) / n * (1-double(iz)/n),
2016 double(iz)/n);
2017 break;
2018 default:
2019 cerr << "case not implementd 878234" << endl;
2020 ploc = 0.0;
2021 }
2022 if (compress[ii] != -1)
2023 locgrid[compress[ii]] = ploc;
2024 }
2025
2026 if (type != TET && type != PRISM) cnt_valid = n3;
2027
2028
2029 if (mesh->GetCurvedElements().IsHighOrder() || 1)
2030 {
2031 mesh->GetCurvedElements().
2032 CalcMultiPointElementTransformation (&locgrid, ei, &grid, &trans);
2033 }
2034 else
2035 {
2036 Vector shape(el.GetNP());
2037 for (int k = 0; k < el.GetNP(); k++)
2038 for (int j = 0; j < 3; j++)
2039 pointmat(k,j) = (*mesh)[el[k]](j);
2040
2041 for (int i = 0; i < cnt_valid; i++)
2042 {
2043 el.GetShapeNew<double> (locgrid[i], shape);
2044 Point<3> pglob;
2045 for (int j = 0; j < 3; j++)
2046 {
2047 pglob(j) = 0;
2048 for (int k = 0; k < el.GetNP(); k++)
2049 pglob(j) += shape(k) * pointmat(k,j);
2050 }
2051 grid[i] = pglob;
2052 }
2053 }
2054
2055 bool has_pos = 0, has_neg = 0;
2056 GetMultiValues( sol, ei, -1, n3,
2057 &locgrid[0](0), &locgrid[1](0)-&locgrid[0](0),
2058 &grid[0](0), &grid[1](0)-&grid[0](0),
2059 &trans[0](0), &trans[1](0)-&trans[0](0),
2060 &val1[0], sol->components);
2061 for (int i = 0; i < cnt_valid; i++)
2062 {
2063 // GetValue (sol, ei, &locgrid[i](0), &grid[i](0), &trans[i](0), comp, val[i]);
2064
2065 // val[i] -= minval;
2066 val1[sol->components*i+comp-1] -= minval;
2067
2068
2069 // if (vsol)
2070 // GetValues (vsol, ei, &locgrid[i](0), &grid[i](0), &trans[i](0), &grads[i](0));
2071 // grads[i] *= -1;
2072
2073 if (val1[i*sol->components+comp-1] > 0)
2074 has_pos = 1;
2075 else
2076 has_neg = 1;
2077 // if (val[i] > 0)
2078 // has_pos = 1;
2079 // else
2080 // has_neg = 1;
2081 }
2082
2083 if (!has_pos || !has_neg) continue;
2084 if (vsol)
2085 {
2086 GetMultiValues(vsol, ei, -1, n3,
2087 &locgrid[0](0), &locgrid[1](0)-&locgrid[0](0),
2088 &grid[0](0), &grid[1](0)-&grid[0](0),
2089 &trans[0](0), &trans[1](0)-&trans[0](0),
2090 &grads1[0](0), vsol->components);
2091 // for (int i = 0; i < cnt_valid; i++)
2092 // grads1[i*sol->components+comp-1] *= -1;
2093 for (int i = 0; i < cnt_valid; i++)
2094 grads1[i] *= -1;
2095
2096 }
2097 for (int ix = 0; ix < n; ix++)
2098 for (int iy = 0; iy < n; iy++)
2099 for (int iz = 0; iz < n; iz++)
2100 {
2101 int base = iz + (n+1)*iy + (n+1)*(n+1)*ix;
2102 int pi[8] =
2103 { base, base+(n+1)*(n+1), base+(n+1)*(n+1)+(n+1), base+(n+1),
2104 base+1, base+(n+1)*(n+1)+1, base+(n+1)*(n+1)+(n+1)+1, base+(n+1)+1 };
2105
2106 for (int j = 0; j < 8; j++)
2107 pi[j] = compress[pi[j]];
2108
2109 int tets[6][4] =
2110 { { 1, 2, 4, 5 },
2111 { 4, 5, 2, 8 },
2112 { 2, 8, 5, 6 },
2113 { 2, 3, 4, 8 },
2114 { 2, 3, 8, 6 },
2115 { 3, 8, 6, 7 } };
2116
2117 for (int ii = 0; ii < 6; ii++)
2118 {
2119 int teti[4];
2120 for (int k = 0; k < 4; k++)
2121 teti[k] = pi[tets[ii][k]-1];
2122
2123 bool is_valid = 1;
2124 for (int j = 0; j < 4; j++)
2125 if (teti[j] == -1) is_valid = 0;
2126
2127 if (!is_valid) continue;
2128
2129 // for (int j = 0; j < 4; j++)
2130 // nodevali[j] = val[teti[j]];
2131 for (int j = 0; j < 4; j++)
2132 nodevali[j] = val1[sol->components*teti[j]+comp-1];
2133
2134 cntce = 0;
2135 for (int j = 0; j < 6; j++)
2136 {
2137 int lpi1 = edgei[j][0];
2138 int lpi2 = edgei[j][1];
2139 if ( (nodevali[lpi1] > 0) !=
2140 (nodevali[lpi2] > 0) )
2141 {
2142 Point<3> p1 = grid[teti[lpi1]];
2143 Point<3> p2 = grid[teti[lpi2]];
2144
2145 edgelam[j] = nodevali[lpi2] / (nodevali[lpi2] - nodevali[lpi1]);
2146 edgep[j] = grid[teti[lpi1]] + (1-edgelam[j]) * (grid[teti[lpi2]]-grid[teti[lpi1]]);
2147 // normp[j] = grads[teti[lpi1]] + (1-edgelam[j]) * (grads[teti[lpi2]]-grads[teti[lpi1]]);
2148 normp[j] = grads1[teti[lpi1]] + (1-edgelam[j]) * (grads1[teti[lpi2]]-grads1[teti[lpi1]]);
2149 // normp[j] = grads1[sol->components*teti[lpi1]+comp-1] + (1-edgelam[j]) * (grads1[sol->components*teti[lpi2]+comp-1]-grads1[sol->components*teti[lpi1]+comp-1]);
2150 cntce++;
2151 cpe3 = cpe2;
2152 cpe2 = cpe1;
2153 cpe1 = j;
2154 if (cntce >= 3)
2155 {
2156 if (!vsol)
2157 {
2158 Point<3> points[3];
2159
2160 points[0] = edgep[cpe1];
2161 points[1] = edgep[cpe2];
2162 points[2] = edgep[cpe3];
2163
2164 Vec<3> normal = Cross (points[2]-points[0], points[1]-points[0]);
2165 if ( ( (normal * (p2-p1)) > 0 ) == ( nodevali[lpi1] < 0) )
2166 normal *= -1;
2167 glNormal3dv (normal);
2168
2169 glVertex3dv (points[0]);
2170 glVertex3dv (points[1]);
2171 glVertex3dv (points[2]);
2172 }
2173 else
2174 {
2175 glNormal3dv (normp[cpe1]);
2176 glVertex3dv (edgep[cpe1]);
2177 glNormal3dv (normp[cpe2]);
2178 glVertex3dv (edgep[cpe2]);
2179 glNormal3dv (normp[cpe3]);
2180 glVertex3dv (edgep[cpe3]);
2181 }
2182 }
2183 }
2184 }
2185 }
2186 }
2187 }
2188 }
2189 glEnd();
2190 }
2191
2192
2193
2194
2195
2196
2197
2198
DrawTrigSurfaceVectors(const NgArray<Point<3>> & lp,const Point<3> & pmin,const Point<3> & pmax,const int sei,const SolData * vsol)2199 void VisualSceneSolution :: DrawTrigSurfaceVectors(const NgArray< Point<3> > & lp,
2200 const Point<3> & pmin, const Point<3> & pmax,
2201 const int sei, const SolData * vsol)
2202 {
2203 shared_ptr<Mesh> mesh = GetMesh();
2204
2205 int dir,dir1,dir2;
2206 double s,t;
2207
2208 Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]);
2209 Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2)));
2210 if (na(0) > na(1) && na(0) > na(2))
2211 dir = 1;
2212 else if (na(1) > na(2))
2213 dir = 2;
2214 else
2215 dir = 3;
2216
2217 dir1 = (dir % 3) + 1;
2218 dir2 = (dir1 % 3) + 1;
2219
2220 Point<2> p2d[3];
2221
2222 int k;
2223
2224 for (k = 0; k < 3; k++)
2225 {
2226 p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad),
2227 (lp[k](dir2-1) - pmin(dir2-1)) / (2*rad));
2228 }
2229
2230
2231 double minx2d, maxx2d, miny2d, maxy2d;
2232 minx2d = maxx2d = p2d[0](0);
2233 miny2d = maxy2d = p2d[0](1);
2234 for (k = 1; k < 3; k++)
2235 {
2236 minx2d = min2 (minx2d, p2d[k](0));
2237 maxx2d = max2 (maxx2d, p2d[k](0));
2238 miny2d = min2 (miny2d, p2d[k](1));
2239 maxy2d = max2 (maxy2d, p2d[k](1));
2240 }
2241
2242 double mat11 = p2d[1](0) - p2d[0](0);
2243 double mat21 = p2d[1](1) - p2d[0](1);
2244 double mat12 = p2d[2](0) - p2d[0](0);
2245 double mat22 = p2d[2](1) - p2d[0](1);
2246
2247 double det = mat11*mat22-mat21*mat12;
2248 double inv11 = mat22/det;
2249 double inv21 = -mat21/det;
2250 double inv12 = -mat12/det;
2251 double inv22 = mat11/det;
2252
2253 // cout << "drawsurfacevectors. xoffset = " << xoffset << ", yoffset = ";
2254 // cout << yoffset << endl;
2255
2256 for (s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize)
2257 if (s >= minx2d && s <= maxx2d)
2258 for (t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize)
2259 if (t >= miny2d && t <= maxy2d)
2260 {
2261 double lam1 = inv11 * (s - p2d[0](0)) + inv12 * (t-p2d[0](1));
2262 double lam2 = inv21 * (s - p2d[0](0)) + inv22 * (t-p2d[0](1));
2263
2264 if (lam1 >= 0 && lam2 >= 0 && lam1+lam2 <= 1)
2265 {
2266 Point<3> cp;
2267 for (k = 0; k < 3; k++)
2268 cp(k) = lp[0](k) +
2269 lam1 * (lp[1](k)-lp[0](k)) +
2270 lam2 * (lp[2](k)-lp[0](k));
2271
2272 Point<2> xref(lam1, lam2);
2273 if (mesh->GetCurvedElements().IsHighOrder())
2274 mesh->GetCurvedElements().
2275 CalcSurfaceTransformation (xref, sei, cp);
2276
2277 Vec<3> v;
2278 double values[6];
2279 bool drawelem =
2280 GetSurfValues (vsol, sei, -1, lam1, lam2, values);
2281
2282 if (!vsol->iscomplex)
2283 for (k = 0; k < 3; k++)
2284 v(k) = values[k];
2285 else
2286 {
2287 if (!imag_part)
2288 for (k = 0; k < 3; k++)
2289 v(k) = values[2*k];
2290 else
2291 for (k = 0; k < 3; k++)
2292 v(k) = values[2*k+1];
2293 }
2294
2295 if (mesh->GetDimension() == 2)
2296 if ( (!vsol->iscomplex && vsol->components != 3) ||
2297 (vsol->iscomplex && vsol->components != 6) )
2298 v(2) = 0;
2299
2300 double val = v.Length();
2301
2302 SetOpenGlColor (val); // (val, minval, maxval, logscale); // change JS
2303
2304 if (val > 1e-10 * maxval)
2305 v *= (rad / val / gridsize * 0.5);
2306 else
2307 drawelem = 0;
2308
2309 if ( drawelem )
2310 DrawCone (cp, cp+4*v, 0.8*rad / gridsize);
2311 }
2312 }
2313
2314 }
2315
2316
2317
DrawSurfaceVectors()2318 void VisualSceneSolution :: DrawSurfaceVectors ()
2319 {
2320 shared_ptr<Mesh> mesh = GetMesh();
2321
2322 SurfaceElementIndex sei;
2323
2324 const SolData * vsol = NULL;
2325 // bool drawelem;
2326
2327 if (vecfunction != -1)
2328 vsol = soldata[vecfunction];
2329
2330 if (mesh->GetTimeStamp () > solutiontimestamp)
2331 vsol = NULL;
2332
2333 if (!vsol) return;
2334
2335
2336 Point<3> pmin = center - Vec3d (rad, rad, rad);
2337 Point<3> pmax = center - Vec3d (rad, rad, rad);
2338
2339
2340 // glColor3d (1.0, 1.0, 1.0);
2341 // glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
2342
2343 if (vsol->draw_surface && showsurfacesolution)
2344 {
2345 int nse = mesh->GetNSE();
2346 for (sei = 0; sei < nse; sei++)
2347 {
2348 const Element2d & el = (*mesh)[sei];
2349
2350 if (el.GetType() == TRIG || el.GetType() == TRIG6)
2351 {
2352
2353 NgArray< Point<3> > lp(3);
2354
2355 lp[0] = mesh->Point(el[2]);
2356 lp[1] = mesh->Point(el[0]);
2357 lp[2] = mesh->Point(el[1]);
2358
2359 DrawTrigSurfaceVectors(lp,pmin,pmax,sei,vsol);
2360
2361 /*
2362 Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]);
2363 Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2)));
2364 if (na(0) > na(1) && na(0) > na(2))
2365 dir = 1;
2366 else if (na(1) > na(2))
2367 dir = 2;
2368 else
2369 dir = 3;
2370
2371 dir1 = (dir % 3) + 1;
2372 dir2 = (dir1 % 3) + 1;
2373
2374 for (k = 0; k < 3; k++)
2375 {
2376 p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad),
2377 (lp[k](dir2-1) - pmin(dir2-1)) / (2*rad));
2378 }
2379
2380 double minx2d, maxx2d, miny2d, maxy2d;
2381 minx2d = maxx2d = p2d[0](0);
2382 miny2d = maxy2d = p2d[0](1);
2383 for (k = 1; k < 3; k++)
2384 {
2385 minx2d = min2 (minx2d, p2d[k](0));
2386 maxx2d = max2 (maxx2d, p2d[k](0));
2387 miny2d = min2 (miny2d, p2d[k](1));
2388 maxy2d = max2 (maxy2d, p2d[k](1));
2389 }
2390
2391 double mat11 = p2d[1](0) - p2d[0](0);
2392 double mat21 = p2d[1](1) - p2d[0](1);
2393 double mat12 = p2d[2](0) - p2d[0](0);
2394 double mat22 = p2d[2](1) - p2d[0](1);
2395
2396 double det = mat11*mat22-mat21*mat12;
2397 double inv11 = mat22/det;
2398 double inv21 = -mat21/det;
2399 double inv12 = -mat12/det;
2400 double inv22 = mat11/det;
2401
2402 // cout << "drawsurfacevectors. xoffset = " << xoffset << ", yoffset = ";
2403 // cout << yoffset << endl;
2404
2405 for (s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize)
2406 if (s >= minx2d && s <= maxx2d)
2407 for (t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize)
2408 if (t >= miny2d && t <= maxy2d)
2409 {
2410 double lam1 = inv11 * (s - p2d[0](0)) + inv12 * (t-p2d[0](1));
2411 double lam2 = inv21 * (s - p2d[0](0)) + inv22 * (t-p2d[0](1));
2412
2413 if (lam1 >= 0 && lam2 >= 0 && lam1+lam2 <= 1)
2414 {
2415 Point<3> cp;
2416 for (k = 0; k < 3; k++)
2417 cp(k) = lp[0](k) +
2418 lam1 * (lp[1](k)-lp[0](k)) +
2419 lam2 * (lp[2](k)-lp[0](k));
2420
2421 Vec<3> v;
2422 double values[6];
2423 drawelem = GetSurfValues (vsol, sei, lam1, lam2, values);
2424
2425 if (!vsol->iscomplex)
2426 for (k = 0; k < 3; k++)
2427 v(k) = values[k];
2428 else
2429 {
2430 if (!imag_part)
2431 for (k = 0; k < 3; k++)
2432 v(k) = values[2*k];
2433 else
2434 for (k = 0; k < 3; k++)
2435 v(k) = values[2*k+1];
2436 }
2437
2438 if (mesh->GetDimension() == 2)
2439 if ( (!vsol->iscomplex && vsol->components != 3) ||
2440 (vsol->iscomplex && vsol->components != 6) )
2441 v(2) = 0;
2442
2443 double val = v.Length();
2444 SetOpenGlColor (val, minval, maxval, logscale);
2445
2446 if (val > 1e-10 * maxval)
2447 v *= (rad / val / gridsize * 0.5);
2448 else drawelem = 0;
2449 // "drawelem": added 07.04.2004 (FB)
2450 if ( drawelem ) DrawCone (cp, cp+4*v, 0.8*rad / gridsize);
2451
2452
2453 }
2454 }
2455 */
2456 }
2457 else if (el.GetType() == QUAD)
2458 {
2459 /*
2460 NgArray < Point<3> > lp(3);
2461
2462 lp[0] = mesh->Point(el[0]);
2463 lp[1] = mesh->Point(el[1]);
2464 lp[2] = mesh->Point(el[2]);
2465
2466 DrawTrigSurfaceVectors(lp,pmin,pmax,sei,vsol);
2467
2468 lp[0] = mesh->Point(el[0]);
2469 lp[1] = mesh->Point(el[2]);
2470 lp[2] = mesh->Point(el[3]);
2471
2472 DrawTrigSurfaceVectors(lp,pmin,pmax,sei,vsol);
2473 */
2474
2475 Point<3> lp[4];
2476 Point<2> p2d[4];
2477
2478 for (int k = 0; k < 4; k++)
2479 lp[k] = mesh->Point (el[k]);
2480
2481
2482 Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]);
2483 Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2)));
2484 int dir, dir1, dir2;
2485 if (na(0) > na(1) && na(0) > na(2))
2486 dir = 1;
2487 else if (na(1) > na(2))
2488 dir = 2;
2489 else
2490 dir = 3;
2491
2492 dir1 = (dir % 3) + 1;
2493 dir2 = (dir1 % 3) + 1;
2494
2495 for (int k = 0; k < 4; k++)
2496 {
2497 p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad),
2498 (lp[k](dir2-1) - pmin(dir2-1)) / (2*rad));
2499 }
2500
2501 double minx2d, maxx2d, miny2d, maxy2d;
2502 minx2d = maxx2d = p2d[0](0);
2503 miny2d = maxy2d = p2d[0](1);
2504 for (int k = 1; k < 4; k++)
2505 {
2506 minx2d = min2 (minx2d, p2d[k](0));
2507 maxx2d = max2 (maxx2d, p2d[k](0));
2508 miny2d = min2 (miny2d, p2d[k](1));
2509 maxy2d = max2 (maxy2d, p2d[k](1));
2510 }
2511
2512 for (double s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize)
2513 if (s >= minx2d && s <= maxx2d)
2514 for (double t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize)
2515 if (t >= miny2d && t <= maxy2d)
2516 {
2517 double lami[3];
2518 Point3d p3d(2*rad*s+pmin(0), 2*rad*t+pmin(1),0);
2519
2520 if (mesh->PointContainedIn2DElement (p3d, lami, sei+1))
2521 {
2522 Point<3> cp = p3d;
2523 double lam1 = lami[0];
2524 double lam2 = lami[1];
2525
2526 //for (k = 0; k < 3; k++)
2527 //cp(k) = lp[0](k) +
2528 //lam1 * (lp[1](k)-lp[0](k)) +
2529 //lam2 * (lp[2](k)-lp[0](k));
2530
2531
2532 Vec<3> v;
2533 double values[6];
2534 bool drawelem = GetSurfValues (vsol, sei, -1, lam1, lam2, values);
2535 (*testout) << "sei " << sei << " lam1 " << lam1 << " lam2 " << lam2 << " drawelem " << drawelem << endl;
2536
2537 if (!vsol->iscomplex)
2538 for (int k = 0; k < 3; k++)
2539 v(k) = values[k];
2540 else
2541 {
2542 if (!imag_part)
2543 for (int k = 0; k < 3; k++)
2544 v(k) = values[2*k];
2545 else
2546 for (int k = 0; k < 3; k++)
2547 v(k) = values[2*k+1];
2548 }
2549
2550 if (mesh->GetDimension() == 2)
2551 if ( (!vsol->iscomplex && vsol->components != 3) ||
2552 (vsol->iscomplex && vsol->components != 6) )
2553 v(2) = 0;
2554
2555 double val = v.Length();
2556 SetOpenGlColor (val); // , minval, maxval, logscale); july 09
2557
2558 (*testout) << "v " << v << endl;
2559
2560 if (val > 1e-10 * maxval)
2561 v *= (rad / val / gridsize * 0.5);
2562
2563 (*testout) << "v " << v << endl;
2564
2565 if ( drawelem )
2566 {
2567 DrawCone (cp, cp+4*v, 0.8*rad / gridsize);
2568 (*testout) << "cp " << cp << " rad " << rad << " gridsize " << gridsize << endl;
2569 }
2570
2571 }
2572 }
2573 }
2574 }
2575 }
2576 }
2577
2578
2579
2580
2581 void VisualSceneSolution ::
DrawIsoLines(const Point<3> & p1,const Point<3> & p2,const Point<3> & p3,double val1,double val2,double val3)2582 DrawIsoLines (const Point<3> & p1,
2583 const Point<3> & p2,
2584 const Point<3> & p3,
2585 double val1, double val2, double val3)
2586 {
2587 DrawIsoLines2 (p1, p2, p1, p3, val1, val2, val1, val3); // , minval, maxval, n);
2588 DrawIsoLines2 (p2, p1, p2, p3, val2, val1, val2, val3); // , minval, maxval, n);
2589 DrawIsoLines2 (p3, p1, p3, p2, val3, val1, val3, val2); // , minval, maxval, n);
2590 }
2591
2592
2593 void VisualSceneSolution ::
DrawIsoLines2(const Point<3> & hp1,const Point<3> & hp2,const Point<3> & hp3,const Point<3> & hp4,double val1,double val2,double val3,double val4)2594 DrawIsoLines2 (const Point<3> & hp1,
2595 const Point<3> & hp2,
2596 const Point<3> & hp3,
2597 const Point<3> & hp4,
2598 double val1, double val2, double val3, double val4)
2599 {
2600 int n = numisolines;
2601 Point<3> p1, p2, p3, p4;
2602 if (val1 < val2)
2603 {
2604 p1 = hp1; p2 = hp2;
2605 }
2606 else
2607 {
2608 p1 = hp2; p2 = hp1;
2609 swap (val1, val2);
2610 }
2611
2612 if (val3 < val4)
2613 {
2614 p3 = hp3; p4 = hp4;
2615 }
2616 else
2617 {
2618 p3 = hp4; p4 = hp3;
2619 swap (val3, val4);
2620 }
2621
2622 val2 += 1e-10;
2623 val4 += 1e-10;
2624
2625 double fac = (maxval-minval) / n;
2626 double idelta1 = 1.0 / (val2 - val1);
2627 double idelta2 = 1.0 / (val4 - val3);
2628
2629 int mini = int ((max2 (val1, val3) - minval) / fac);
2630 int maxi = int ((min2 (val2, val4) - minval) / fac);
2631 if (mini < 0) mini = 0;
2632 if (maxi > n-1) maxi = n-1;
2633
2634 for (int i = mini; i <= maxi; i++)
2635 {
2636 double val = minval + i * fac;
2637 double lam1 = (val - val1) * idelta1;
2638 double lam2 = (val - val3) * idelta2;
2639 if (lam1 >= 0 && lam1 <= 1 && lam2 >= 0 && lam2 <= 1)
2640 {
2641 Point<3> lp1 = p1 + lam1 * (p2-p1);
2642 Point<3> lp2 = p3 + lam2 * (p4-p3);
2643 glVertex3dv (lp1 );
2644 glVertex3dv (lp2 );
2645 // glVertex3dv (lp2 ); // better ?
2646 // glVertex3dv (lp1 );
2647 }
2648 }
2649 }
2650
2651
2652
2653 void VisualSceneSolution ::
GetMinMax(int funcnr,int comp,double & minv,double & maxv) const2654 GetMinMax (int funcnr, int comp, double & minv, double & maxv) const
2655 {
2656 shared_ptr<Mesh> mesh = GetMesh();
2657
2658 // static int timer1 = NgProfiler::CreateTimer ("getminmax, vol");
2659 // static int timer2 = NgProfiler::CreateTimer ("getminmax, surf");
2660
2661 #ifdef PARALLEL
2662 auto comm = mesh->GetCommunicator();
2663 if (comm.Size() > 1)
2664 {
2665 if (id == 0)
2666 {
2667 MyMPI_SendCmd ("redraw");
2668 MyMPI_SendCmd ("getminmax");
2669 }
2670 MyMPI_Bcast (funcnr, mesh->GetCommunicator());
2671 MyMPI_Bcast (comp, mesh->GetCommunicator());
2672 }
2673 #endif
2674
2675 // double val;
2676 // bool considerElem;
2677
2678 bool hasit = false;
2679 #ifdef max
2680 #undef max
2681 #endif
2682 minv = numeric_limits<double>::max();
2683 maxv = -numeric_limits<double>::max();
2684
2685 if ((ntasks == 1) || (id > 0))
2686 if (funcnr != -1)
2687 {
2688 const SolData * sol = soldata[funcnr];
2689
2690 if (sol->draw_volume)
2691 {
2692 // NgProfiler::RegionTimer reg1 (timer1);
2693
2694 int ne = mesh->GetNE();
2695
2696 mutex min_mutex;
2697 mutex max_mutex;
2698
2699 ParallelFor(0, ne, [&] (int first, int next)
2700 {
2701 double minv_local = numeric_limits<double>::max();
2702 double maxv_local = -numeric_limits<double>::max();
2703 for (int i=first; i<next; i++)
2704 {
2705 double val;
2706 bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val);
2707 if (considerElem)
2708 {
2709 if (val > maxv_local) maxv_local = val;
2710 if (val < minv_local) minv_local = val;
2711 hasit = true;
2712 }
2713 }
2714 if(minv_local < minv)
2715 {
2716 lock_guard<mutex> guard(min_mutex);
2717 if(minv_local < minv)
2718 minv = minv_local;
2719 }
2720 if(maxv_local > maxv)
2721 {
2722 lock_guard<mutex> guard(max_mutex);
2723 if(maxv_local > maxv)
2724 maxv = maxv_local;
2725 }
2726 });
2727 }
2728
2729 if (sol->draw_surface)
2730 {
2731 // NgProfiler::RegionTimer reg2 (timer2);
2732
2733 // int nse = mesh->GetNSE();
2734 // for (int i = 0; i < nse; i++)
2735 for (SurfaceElementIndex i : mesh->SurfaceElements().Range())
2736 {
2737 ELEMENT_TYPE type = (*mesh)[i].GetType();
2738 double val;
2739 bool considerElem = (type == QUAD)
2740 ? GetSurfValue (sol, i, -1, 0.5, 0.5, comp, val)
2741 : GetSurfValue (sol, i, -1, 0.3333333, 0.3333333, comp, val);
2742
2743 if (considerElem)
2744 {
2745 if (val > maxv) maxv = val;
2746 if (val < minv) minv = val;
2747 hasit = true;
2748 }
2749 }
2750 }
2751 }
2752 if (minv == maxv) maxv = minv+1e-6;
2753 if (!hasit) { minv = 0; maxv = 1; }
2754
2755 #ifdef PARALLEL
2756 if ((ntasks > 1) && (id == 0))
2757 {
2758 minv = 1e99;
2759 maxv = -1e99;
2760 }
2761 if (ntasks > 1)
2762 {
2763 double hmin, hmax;
2764 MPI_Reduce (&minv, &hmin, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
2765 MPI_Reduce (&maxv, &hmax, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
2766 minv = hmin;
2767 maxv = hmax;
2768 }
2769 #endif
2770 }
2771
2772
2773
2774
2775
2776 bool VisualSceneSolution ::
GetValues(const SolData * data,ElementIndex elnr,double lam1,double lam2,double lam3,double * values) const2777 GetValues (const SolData * data, ElementIndex elnr,
2778 double lam1, double lam2, double lam3,
2779 double * values) const
2780 {
2781 bool ok = false;
2782 switch (data->soltype)
2783 {
2784 case SOL_VIRTUALFUNCTION:
2785 {
2786 ok = data->solclass->GetValue (elnr, lam1, lam2, lam3, values);
2787 break;
2788 }
2789 default:
2790 {
2791 for (int i = 0; i < data->components; i++)
2792 ok = GetValue (data, elnr, lam1, lam2, lam3, i+1, values[i]);
2793 }
2794 }
2795 return ok;
2796 }
2797
2798 bool VisualSceneSolution ::
GetValues(const SolData * data,ElementIndex elnr,const double xref[],const double x[],const double dxdxref[],double * values) const2799 GetValues (const SolData * data, ElementIndex elnr,
2800 const double xref[], const double x[], const double dxdxref[],
2801 double * values) const
2802 {
2803 bool ok = false;
2804 switch (data->soltype)
2805 {
2806 case SOL_VIRTUALFUNCTION:
2807 {
2808 ok = data->solclass->GetValue (elnr, xref, x, dxdxref, values);
2809 break;
2810 }
2811 default:
2812 {
2813 for (int i = 0; i < data->components; i++)
2814 ok = GetValue (data, elnr, xref[0], xref[1], xref[2], i+1, values[i]);
2815 }
2816 }
2817 return ok;
2818 }
2819
2820
2821 bool VisualSceneSolution ::
GetValue(const SolData * data,ElementIndex elnr,const double xref[],const double x[],const double dxdxref[],int comp,double & val) const2822 GetValue (const SolData * data, ElementIndex elnr,
2823 const double xref[], const double x[], const double dxdxref[],
2824 int comp, double & val) const
2825 {
2826 shared_ptr<Mesh> mesh = GetMesh();
2827
2828 double lam1 = xref[0];
2829 double lam2 = xref[1];
2830 double lam3 = xref[2];
2831
2832 val = 0;
2833 bool ok = 0;
2834
2835
2836 if (comp == 0)
2837 {
2838 NgArrayMem<double,20> values(data->components);
2839 ok = GetValues (data, elnr, xref, x, dxdxref, &values[0]);
2840
2841 val = ExtractValue (data, 0, &values[0]);
2842 return ok;
2843 }
2844
2845
2846 switch (data->soltype)
2847 {
2848 case SOL_VIRTUALFUNCTION:
2849 {
2850 double values[20];
2851 ok = data->solclass->GetValue (elnr, xref, x, dxdxref, values);
2852
2853 val = values[comp-1];
2854 return ok;
2855 }
2856 case SOL_NODAL:
2857 {
2858 const Element & el = (*mesh)[elnr];
2859
2860 double lami[8] = { 0.0 };
2861 int np = 0;
2862
2863 switch (el.GetType())
2864 {
2865 case TET:
2866 case TET10:
2867 {
2868 lami[1] = lam1;
2869 lami[2] = lam2;
2870 lami[3] = lam3;
2871 lami[0] = 1-lam1-lam2-lam3;
2872 np = 4;
2873 break;
2874 }
2875 case PRISM:
2876 case PRISM12:
2877 case PRISM15:
2878 {
2879 lami[0] = (1-lam3) * (1-lam1-lam2);
2880 lami[1] = (1-lam3) * lam1;
2881 lami[2] = (1-lam3) * lam2;
2882 lami[3] = (lam3) * (1-lam1-lam2);
2883 lami[4] = (lam3) * lam1;
2884 lami[5] = (lam3) * lam2;
2885 np = 6;
2886 break;
2887 }
2888 default:
2889 cerr << "case not implementd 23523" << endl;
2890 }
2891
2892 for (int i = 0; i < np; i++)
2893 val += lami[i] * data->data[(el[i]-1) * data->dist + comp-1];
2894
2895 return 1;
2896 }
2897
2898 case SOL_ELEMENT:
2899 {
2900 val = data->data[elnr * data->dist + comp-1];
2901 return 1;
2902 }
2903
2904 case SOL_SURFACE_ELEMENT:
2905 return 0;
2906
2907 case SOL_NONCONTINUOUS:
2908 {
2909 const Element & el = (*mesh)[elnr];
2910
2911 double lami[8] = { 0.0 };
2912 int np = 0;
2913
2914 switch (el.GetType())
2915 {
2916 case TET:
2917 case TET10:
2918 {
2919 lami[1] = lam1;
2920 lami[2] = lam2;
2921 lami[3] = lam3;
2922 lami[0] = 1-lam1-lam2-lam3;
2923 np = 4;
2924 break;
2925 }
2926 case PRISM:
2927 case PRISM12:
2928 case PRISM15:
2929 {
2930 lami[0] = (1-lam3) * (1-lam1-lam2);
2931 lami[1] = (1-lam3) * lam1;
2932 lami[2] = (1-lam3) * lam2;
2933 lami[3] = (lam3) * (1-lam1-lam2);
2934 lami[4] = (lam3) * lam1;
2935 lami[5] = (lam3) * lam2;
2936 np = 6;
2937 break;
2938 }
2939 case PYRAMID:
2940 case PYRAMID13:
2941 {
2942 if (lam3 > 1-1e-5)
2943 {
2944 lami[0] = lami[1] = lami[2] = lami[3] = 0;
2945 lami[4] = 1;
2946 }
2947 else
2948 {
2949 double x0 = lam1 / (1-lam3);
2950 double y0 = lam2 / (1-lam3);
2951 lami[0] = (1-x0) * (1-y0) * (1-lam3);
2952 lami[1] = ( x0) * (1-y0) * (1-lam3);
2953 lami[2] = ( x0) * ( y0) * (1-lam3);
2954 lami[3] = (1-x0) * ( y0) * (1-lam3);
2955 lami[4] = lam3;
2956 np = 5;
2957 }
2958 break;
2959 }
2960 default:
2961 np = 0;
2962 }
2963
2964 int base;
2965 if (data->order == 1)
2966 base = 6 * elnr;
2967 else
2968 base = 10 * elnr;
2969
2970
2971 for (int i = 0; i < np; i++)
2972 val += lami[i] * data->data[(base+i) * data->dist + comp-1];
2973
2974 return 1;
2975 }
2976
2977 case SOL_MARKED_ELEMENTS:
2978 {
2979 val = (*mesh)[elnr].TestRefinementFlag();
2980 return 1;
2981 }
2982
2983 case SOL_ELEMENT_ORDER:
2984 {
2985 val = (*mesh)[elnr].GetOrder();
2986 return 1;
2987 }
2988
2989 default:
2990 cerr << "case not handled 7234" << endl;
2991 }
2992 return 0;
2993 }
2994
2995
2996
2997 bool VisualSceneSolution ::
GetValue(const SolData * data,ElementIndex elnr,double lam1,double lam2,double lam3,int comp,double & val) const2998 GetValue (const SolData * data, ElementIndex elnr,
2999 double lam1, double lam2, double lam3,
3000 int comp, double & val) const
3001 {
3002 shared_ptr<Mesh> mesh = GetMesh();
3003
3004 val = 0;
3005 bool ok = 0;
3006
3007 if (comp == 0)
3008 {
3009 NgArrayMem<double,20> values(data->components);
3010 ok = GetValues (data, elnr, lam1, lam2, lam3, &values[0]);
3011 val = ExtractValue (data, 0, &values[0]);
3012 return ok;
3013 }
3014
3015
3016 switch (data->soltype)
3017 {
3018 case SOL_VIRTUALFUNCTION:
3019 {
3020 val = 0.0;
3021 double values[20];
3022 ok = data->solclass->GetValue (elnr, lam1, lam2, lam3, values);
3023
3024 val = values[comp-1];
3025 return ok;
3026 }
3027 case SOL_NODAL:
3028 {
3029 const Element & el = (*mesh)[elnr];
3030
3031 double lami[8] = { 0.0 };
3032 int np = 0;
3033
3034 switch (el.GetType())
3035 {
3036 case TET:
3037 case TET10:
3038 {
3039 lami[1] = lam1;
3040 lami[2] = lam2;
3041 lami[3] = lam3;
3042 lami[0] = 1-lam1-lam2-lam3;
3043 np = 4;
3044 break;
3045 }
3046 case PRISM:
3047 case PRISM12:
3048 {
3049 lami[0] = (1-lam3) * (1-lam1-lam2);
3050 lami[1] = (1-lam3) * lam1;
3051 lami[2] = (1-lam3) * lam2;
3052 lami[3] = (lam3) * (1-lam1-lam2);
3053 lami[4] = (lam3) * lam1;
3054 lami[5] = (lam3) * lam2;
3055 np = 6;
3056 break;
3057 }
3058 default:
3059 cerr << "case not implemented 234324" << endl;
3060 }
3061
3062 for (int i = 0; i < np; i++)
3063 val += lami[i] * data->data[(el[i]-1) * data->dist + comp-1];
3064
3065 return 1;
3066 }
3067
3068 case SOL_ELEMENT:
3069 {
3070 val = data->data[elnr * data->dist + comp-1];
3071 return 1;
3072 }
3073
3074 case SOL_SURFACE_ELEMENT:
3075 return 0;
3076
3077 case SOL_NONCONTINUOUS:
3078 {
3079 const Element & el = (*mesh)[elnr];
3080
3081 double lami[8] = { 0.0 };
3082 int np = 0;
3083
3084 switch (el.GetType())
3085 {
3086 case TET:
3087 case TET10:
3088 {
3089 lami[1] = lam1;
3090 lami[2] = lam2;
3091 lami[3] = lam3;
3092 lami[0] = 1-lam1-lam2-lam3;
3093 np = 4;
3094 break;
3095 }
3096 case PRISM:
3097 case PRISM12:
3098 {
3099 lami[0] = (1-lam3) * (1-lam1-lam2);
3100 lami[1] = (1-lam3) * lam1;
3101 lami[2] = (1-lam3) * lam2;
3102 lami[3] = (lam3) * (1-lam1-lam2);
3103 lami[4] = (lam3) * lam1;
3104 lami[5] = (lam3) * lam2;
3105 np = 6;
3106 break;
3107 }
3108 case PYRAMID:
3109 {
3110 if (lam3 > 1-1e-5)
3111 {
3112 lami[0] = lami[1] = lami[2] = lami[3] = 0;
3113 lami[4] = 1;
3114 }
3115 else
3116 {
3117 double x0 = lam1 / (1-lam3);
3118 double y0 = lam2 / (1-lam3);
3119 lami[0] = (1-x0) * (1-y0) * (1-lam3);
3120 lami[1] = ( x0) * (1-y0) * (1-lam3);
3121 lami[2] = ( x0) * ( y0) * (1-lam3);
3122 lami[3] = (1-x0) * ( y0) * (1-lam3);
3123 lami[4] = lam3;
3124 np = 5;
3125 }
3126 break;
3127 }
3128 default:
3129 np = 0;
3130 }
3131
3132 int base;
3133 if (data->order == 1)
3134 base = 6 * elnr;
3135 else
3136 base = 10 * elnr;
3137
3138
3139 for (int i = 0; i < np; i++)
3140 val += lami[i] * data->data[(base+i) * data->dist + comp-1];
3141
3142 return 1;
3143 }
3144
3145 case SOL_MARKED_ELEMENTS:
3146 {
3147 val = (*mesh)[elnr].TestRefinementFlag();
3148 return 1;
3149 }
3150
3151 case SOL_ELEMENT_ORDER:
3152 {
3153 val = (*mesh)[elnr].GetOrder();
3154 return 1;
3155 }
3156 default:
3157 cerr << "case not implemented 234234" << endl;
3158 }
3159 return 0;
3160 }
3161
3162
3163
3164
3165
3166
3167
3168 bool VisualSceneSolution ::
GetValueComplex(const SolData * data,ElementIndex elnr,double lam1,double lam2,double lam3,int comp,complex<double> & val) const3169 GetValueComplex (const SolData * data, ElementIndex elnr,
3170 double lam1, double lam2, double lam3,
3171 int comp, complex<double> & val) const
3172 {
3173 shared_ptr<Mesh> mesh = GetMesh();
3174
3175 val = 0.0;
3176 bool ok = 0;
3177
3178
3179 switch (data->soltype)
3180 {
3181 case SOL_VIRTUALFUNCTION:
3182 {
3183 double values[20];
3184 ok = data->solclass->GetValue (elnr, lam1, lam2, lam3, values);
3185 val = complex<double> (values[comp-1], values[comp]);
3186 return ok;
3187 }
3188 default:
3189 cerr << "case not handled 234234" << endl;
3190 }
3191 return 0;
3192 }
3193
3194
3195 bool VisualSceneSolution ::
GetMultiValues(const SolData * data,ElementIndex elnr,int facetnr,int npt,const double * xref,int sxref,const double * x,int sx,const double * dxdxref,int sdxdxref,double * val,int sval) const3196 GetMultiValues (const SolData * data, ElementIndex elnr, int facetnr, int npt,
3197 const double * xref, int sxref,
3198 const double * x, int sx,
3199 const double * dxdxref, int sdxdxref,
3200 double * val, int sval) const
3201 {
3202 bool drawelem = false;
3203 if (data->soltype == SOL_VIRTUALFUNCTION)
3204 drawelem = data->solclass->GetMultiValue(elnr, facetnr, npt, xref, sxref, x, sx, dxdxref, sdxdxref, val, sval);
3205 else
3206 for (int i = 0; i < npt; i++)
3207 drawelem = GetValues (data, elnr, xref+i*sxref, x+i*sx, dxdxref+i*sdxdxref, val+i*sval);
3208 return drawelem;
3209 }
3210
3211
3212
3213
3214
3215
3216 bool VisualSceneSolution ::
GetSurfValues(const SolData * data,SurfaceElementIndex selnr,int facetnr,double lam1,double lam2,double * values) const3217 GetSurfValues (const SolData * data, SurfaceElementIndex selnr, int facetnr,
3218 double lam1, double lam2,
3219 double * values) const
3220 {
3221 bool ok = false;
3222 switch (data->soltype)
3223 {
3224 case SOL_VIRTUALFUNCTION:
3225 {
3226 ok = data->solclass->GetSurfValue (selnr, facetnr, lam1, lam2, values);
3227 // ok = 1;
3228 // values[0] = 1.0;
3229 break;
3230 }
3231 default:
3232 {
3233 for (int i = 0; i < data->components; i++)
3234 ok = GetSurfValue (data, selnr, facetnr, lam1, lam2, i+1, values[i]);
3235 }
3236 }
3237 return ok;
3238 }
3239
3240
3241 bool VisualSceneSolution ::
GetSurfValues(const SolData * data,SurfaceElementIndex selnr,int facetnr,const double xref[],const double x[],const double dxdxref[],double * values) const3242 GetSurfValues (const SolData * data, SurfaceElementIndex selnr, int facetnr,
3243 const double xref[], const double x[], const double dxdxref[],
3244 double * values) const
3245 {
3246 bool ok = false;
3247 switch (data->soltype)
3248 {
3249 case SOL_VIRTUALFUNCTION:
3250 {
3251 ok = data->solclass->GetSurfValue (selnr, facetnr, xref, x, dxdxref, values);
3252 break;
3253 }
3254 default:
3255 {
3256 for (int i = 0; i < data->components; i++)
3257 ok = GetSurfValue (data, selnr, facetnr, xref[0], xref[1], i+1, values[i]);
3258 }
3259 }
3260 return ok;
3261 }
3262
3263 bool VisualSceneSolution ::
GetMultiSurfValues(const SolData * data,SurfaceElementIndex elnr,int facetnr,int npt,const double * xref,int sxref,const double * x,int sx,const double * dxdxref,int sdxdxref,double * val,int sval) const3264 GetMultiSurfValues (const SolData * data, SurfaceElementIndex elnr, int facetnr, int npt,
3265 const double * xref, int sxref,
3266 const double * x, int sx,
3267 const double * dxdxref, int sdxdxref,
3268 double * val, int sval) const
3269 {
3270 bool drawelem = false;
3271 if (data->soltype == SOL_VIRTUALFUNCTION)
3272 drawelem = data->solclass->GetMultiSurfValue(elnr, facetnr, npt, xref, sxref, x, sx, dxdxref, sdxdxref, val, sval);
3273 else
3274 for (int i = 0; i < npt; i++)
3275 drawelem = GetSurfValues (data, elnr, facetnr, xref+i*sxref, x+i*sx, dxdxref+i*sdxdxref, val+i*sval);
3276 return drawelem;
3277 }
3278
ExtractValue(const SolData * data,int comp,double * values) const3279 double VisualSceneSolution :: ExtractValue (const SolData * data, int comp, double * values) const
3280 {
3281 double val = 0;
3282 if (comp == 0)
3283 {
3284 switch (evalfunc)
3285 {
3286 case FUNC_ABS:
3287 {
3288 for (int ci = 0; ci < data->components; ci++)
3289 val += sqr (values[ci]);
3290 val = sqrt (val);
3291 break;
3292 }
3293 case FUNC_ABS_TENSOR:
3294 {
3295 int d = 0;
3296 switch (data->components)
3297 {
3298 case 1: d = 1; break;
3299 case 3: d = 2; break;
3300 case 6: d = 3; break;
3301 }
3302 for (int ci = 0; ci < d; ci++)
3303 val += sqr (values[ci]);
3304 for (int ci = d; ci < data->components; ci++)
3305 val += 2*sqr (values[ci]);
3306 val = sqrt (val);
3307 break;
3308 }
3309
3310 case FUNC_MISES:
3311 {
3312 int d = 0;
3313 switch(data->components)
3314 {
3315 case 1: d = 1; break;
3316 case 3: d = 2; break;
3317 case 6: d = 3; break;
3318 }
3319 int ci;
3320 double trace = 0.;
3321 for (ci = 0; ci < d; ci++)
3322 trace += 1./3.*(values[ci]);
3323 for (ci = 0; ci < d; ci++)
3324 val += sqr (values[ci]-trace);
3325 for (ci = d; ci < data->components; ci++)
3326 val += 2.*sqr (values[ci]);
3327 val = sqrt (val);
3328 break;
3329 }
3330 case FUNC_MAIN:
3331 {
3332 int d = 0;
3333 switch(data->components)
3334 {
3335 case 1: d = 1; break;
3336 case 3: d = 2; break;
3337 case 6: d = 3; break;
3338 }
3339 Mat<3,3> m ;
3340 Vec<3> ev;
3341 int ci;
3342 for (ci = 0; ci < d; ci++)
3343 m(ci,ci) = (values[ci]);
3344 m(0,1) = m(1,0) = values[3];
3345 m(0,2) = m(2,0) = values[4];
3346 m(1,2) = m(2,1) = values[5];
3347
3348 EigenValues (m, ev);
3349 double help;
3350 for (int i=0; i<d; i++)
3351 {
3352 for (int j=d-1; i<j; j--)
3353 {
3354 if ( abs(ev(j)) > abs(ev(j-1)) )
3355 {
3356 help = ev(j);
3357 ev(j) = ev(j-1);
3358 ev(j-1) = help;
3359 }
3360 }
3361 }
3362 val = (ev(0));
3363 break;
3364 }
3365 }
3366 return val;
3367 }
3368
3369 return values[comp-1];
3370 }
3371
ExtractValueComplex(const SolData * data,int comp,double * values) const3372 complex<double> VisualSceneSolution :: ExtractValueComplex (const SolData * data, int comp, double * values) const
3373 {
3374 if (!data->iscomplex)
3375 return values[comp-1];
3376 else
3377 return complex<double> (values[comp-1], values[comp]);
3378 }
3379
3380
3381
3382
3383 bool VisualSceneSolution ::
GetSurfValueComplex(const SolData * data,SurfaceElementIndex selnr,int facetnr,double lam1,double lam2,int comp,complex<double> & val) const3384 GetSurfValueComplex (const SolData * data, SurfaceElementIndex selnr, int facetnr,
3385 double lam1, double lam2,
3386 int comp, complex<double> & val) const
3387 {
3388 switch (data->soltype)
3389 {
3390 case SOL_VIRTUALFUNCTION:
3391 {
3392 NgArrayMem<double,20> values(data->components);
3393 bool ok;
3394
3395 ok = data->solclass->GetSurfValue (selnr, facetnr, lam1, lam2, &values[0]);
3396
3397 if (ok)
3398 {
3399 if (!data->iscomplex)
3400 val = values[comp-1];
3401 else
3402 val = complex<double> (values[comp-1], values[comp]);
3403 }
3404
3405 return ok;
3406 }
3407 default:
3408 cerr << "case not implementd 6565" << endl;
3409 }
3410 return 0;
3411 }
3412
3413 bool VisualSceneSolution ::
GetSurfValue(const SolData * data,SurfaceElementIndex selnr,int facetnr,double lam1,double lam2,int comp,double & val) const3414 GetSurfValue (const SolData * data, SurfaceElementIndex selnr, int facetnr,
3415 double lam1, double lam2,
3416 int comp, double & val) const
3417 {
3418 bool ok;
3419 if (comp == 0)
3420 {
3421 val = 0;
3422 NgArrayMem<double,20> values(data->components);
3423 ok = GetSurfValues (data, selnr, facetnr, lam1, lam2, &values[0]);
3424 val = ExtractValue (data, 0, &values[0]);
3425 return ok;
3426 }
3427
3428
3429 switch (data->soltype)
3430 {
3431 case SOL_VIRTUALFUNCTION:
3432 {
3433
3434 NgArrayMem<double,20> values(data->components);
3435 bool ok;
3436
3437 ok = data->solclass->GetSurfValue (selnr, facetnr, lam1, lam2, &values[0]);
3438
3439 if (ok)
3440 {
3441 if (!data->iscomplex)
3442 val = values[comp-1];
3443 else
3444 {
3445 // cout << "time = " << time << ", cos = " << cos(time) << endl;
3446
3447 // old version: val = values[comp-1]*cos(3*time) + values[comp]*sin(3*time);
3448 // SZ: Sept 06
3449 if(comp%2==0)
3450 val = values[comp-1]*cos(3*time) - values[comp-2]*sin(3*time);
3451 else
3452 val = values[comp-1]*cos(3*time) + values[comp]*sin(3*time);
3453
3454
3455
3456 }
3457 }
3458
3459 return ok;
3460 }
3461
3462
3463 case SOL_NODAL:
3464 {
3465 shared_ptr<Mesh> mesh = GetMesh();
3466 const Element2d & el = (*mesh)[selnr];
3467
3468 double lami[8];
3469 int np, i;
3470 val = 0;
3471 double lam3 = 1-lam1-lam2;
3472
3473 switch (el.GetType())
3474 {
3475 case TRIG:
3476 /*
3477 lami[0] = lam3;
3478 lami[1] = lam1;
3479 lami[2] = lam2;
3480 */
3481 lami[0] = lam1;
3482 lami[1] = lam2;
3483 lami[2] = lam3;
3484 np = 3;
3485 break;
3486
3487 case TRIG6:
3488 /*
3489 lami[0] = lam3*(2*lam3-1);
3490 lami[1] = lam1*(2*lam1-1);
3491 lami[2] = lam2*(2*lam2-1);
3492 */
3493 // hierarchical basis:
3494 lami[0] = lam3;
3495 lami[1] = lam1;
3496 lami[2] = lam2;
3497 lami[3] = 4*lam1*lam2;
3498 lami[4] = 4*lam2*lam3;
3499 lami[5] = 4*lam1*lam3;
3500 np = 6;
3501 break;
3502
3503 case QUAD:
3504 case QUAD6:
3505 case QUAD8:
3506 lami[0] = (1-lam1)*(1-lam2);
3507 lami[1] = lam1 * (1-lam2);
3508 lami[2] = lam1 * lam2;
3509 lami[3] = (1-lam1) * lam2;
3510 np = 4;
3511 break;
3512
3513 default:
3514 np = 0;
3515 }
3516
3517 for (i = 0; i < np; i++)
3518 val += lami[i] * data->data[(el[i]-1) * data->dist + comp-1];
3519
3520 return 1;
3521 }
3522
3523 case SOL_ELEMENT:
3524 {
3525 shared_ptr<Mesh> mesh = GetMesh();
3526 int el1, el2;
3527 mesh->GetTopology().GetSurface2VolumeElement (selnr+1, el1, el2);
3528 el1--;
3529
3530 val = data->data[el1 * data->dist+comp-1];
3531 return 1;
3532 }
3533
3534 case SOL_NONCONTINUOUS:
3535 {
3536 val = 0;
3537 // ?????
3538 return 0;
3539 }
3540
3541 case SOL_SURFACE_ELEMENT:
3542 {
3543 val = data->data[selnr * data->dist + comp-1];
3544 return 1;
3545 }
3546
3547 case SOL_SURFACE_NONCONTINUOUS:
3548 {
3549 shared_ptr<Mesh> mesh = GetMesh();
3550 const Element2d & el = (*mesh)[selnr];
3551
3552 double lami[8];
3553 int np = 0;
3554 val = 0;
3555 int order = data->order;
3556
3557 switch (order)
3558 {
3559 case 0:
3560 return data->data[selnr * data->dist + comp-1];
3561 case 1:
3562 {
3563 switch (el.GetType())
3564 {
3565 case TRIG:
3566 case TRIG6:
3567 {
3568 lami[1] = lam1;
3569 lami[2] = lam2;
3570 lami[0] = 1-lam1-lam2;
3571 np = 3;
3572 break;
3573 }
3574 default:
3575 cerr << "case not implementd 2342" << endl;
3576 }
3577 break;
3578 }
3579 case 2:
3580 {
3581 switch (el.GetType())
3582 {
3583 case TRIG:
3584 {
3585 lami[1] = lam1;
3586 lami[2] = lam2;
3587 lami[0] = 1-lam1-lam2;
3588 np = 3;
3589 break;
3590 }
3591 case TRIG6:
3592 {
3593 double lam3 = 1-lam1-lam2;
3594 lami[1] = 2*lam1 * (lam1-0.5);
3595 lami[2] = 2*lam2 * (lam2-0.5);
3596 lami[0] = 2*lam3 * (lam3-0.5);
3597 lami[3] = 4*lam1*lam2;
3598 lami[4] = 4*lam2*lam3;
3599 lami[5] = 4*lam1*lam3;
3600 np = 6;
3601 break;
3602 }
3603 default:
3604 cerr << "case not implemented 8712" << endl;
3605 }
3606 break;
3607 }
3608 }
3609
3610 int base;
3611 if (order == 1)
3612 base = 4 * selnr;
3613 else
3614 base = 9 * selnr;
3615
3616 for (int i = 0; i < np; i++)
3617 val += lami[i] * data->data[(base+i) * data->dist + comp-1];
3618
3619 return 1;
3620 }
3621
3622 case SOL_MARKED_ELEMENTS:
3623 {
3624 shared_ptr<Mesh> mesh = GetMesh();
3625 val = (*mesh)[selnr].TestRefinementFlag();
3626 return 1;
3627 }
3628
3629 case SOL_ELEMENT_ORDER:
3630 {
3631 shared_ptr<Mesh> mesh = GetMesh();
3632 val = (*mesh)[selnr].GetOrder();
3633 return 1;
3634 }
3635
3636 }
3637 return 0;
3638 }
3639
3640
3641
3642
3643
3644
3645
3646
3647
3648
3649
3650
3651 bool VisualSceneSolution ::
GetSurfValue(const SolData * data,SurfaceElementIndex selnr,int facetnr,const double xref[],const double x[],const double dxdxref[],int comp,double & val) const3652 GetSurfValue (const SolData * data, SurfaceElementIndex selnr, int facetnr,
3653 const double xref[], const double x[], const double dxdxref[],
3654 int comp, double & val) const
3655 {
3656 shared_ptr<Mesh> mesh = GetMesh();
3657
3658 double lam1 = xref[0], lam2 = xref[1];
3659
3660 bool ok;
3661 if (comp == 0)
3662 {
3663 val = 0;
3664 NgArrayMem<double,20> values(data->components);
3665 ok = GetSurfValues (data, selnr, facetnr, xref, x, dxdxref, &values[0]);
3666 val = ExtractValue (data, 0, &values[0]);
3667 return ok;
3668 }
3669
3670
3671 switch (data->soltype)
3672 {
3673 case SOL_VIRTUALFUNCTION:
3674 {
3675 NgArrayMem<double,20> values(data->components);
3676 bool ok;
3677
3678 // ok = data->solclass->GetSurfValue (selnr, lam1, lam2, &values[0]);
3679 // cout << "data->solclass = " << flush << data->solclass << endl;
3680 ok = data->solclass->GetSurfValue (selnr, facetnr, xref, x, dxdxref, &values[0]);
3681 // ok = 1;
3682 // values[0] = 1.0;
3683
3684 if (ok)
3685 {
3686 if (!data->iscomplex)
3687 val = values[comp-1];
3688 else
3689 {
3690 // cout << "time = " << time << ", cos = " << cos(time) << endl;
3691
3692 // old version: val = values[comp-1]*cos(3*time) + values[comp]*sin(3*time);
3693 // SZ: Sept 06
3694 if(comp%2==0)
3695 val = values[comp-1]*cos(3*time) - values[comp-2]*sin(3*time);
3696 else
3697 val = values[comp-1]*cos(3*time) + values[comp]*sin(3*time);
3698
3699 }
3700 }
3701
3702 return ok;
3703 }
3704
3705
3706 case SOL_NODAL:
3707 {
3708 const Element2d & el = (*mesh)[selnr];
3709
3710 double lami[8];
3711 int np, i;
3712 val = 0;
3713 double lam3 = 1-lam1-lam2;
3714
3715 switch (el.GetType())
3716 {
3717 case TRIG:
3718 /*
3719 lami[0] = lam3;
3720 lami[1] = lam1;
3721 lami[2] = lam2;
3722 */
3723 lami[0] = lam1;
3724 lami[1] = lam2;
3725 lami[2] = lam3;
3726 np = 3;
3727 break;
3728
3729 case TRIG6:
3730 /*
3731 lami[0] = lam3*(2*lam3-1);
3732 lami[1] = lam1*(2*lam1-1);
3733 lami[2] = lam2*(2*lam2-1);
3734 */
3735 // hierarchical basis:
3736 lami[0] = lam3;
3737 lami[1] = lam1;
3738 lami[2] = lam2;
3739 lami[3] = 4*lam1*lam2;
3740 lami[4] = 4*lam2*lam3;
3741 lami[5] = 4*lam1*lam3;
3742 np = 6;
3743 break;
3744
3745 case QUAD:
3746 case QUAD6:
3747 case QUAD8:
3748 lami[0] = (1-lam1)*(1-lam2);
3749 lami[1] = lam1 * (1-lam2);
3750 lami[2] = lam1 * lam2;
3751 lami[3] = (1-lam1) * lam2;
3752 np = 4;
3753 break;
3754
3755 default:
3756 np = 0;
3757 }
3758
3759 for (i = 0; i < np; i++)
3760 val += lami[i] * data->data[(el[i]-1) * data->dist + comp-1];
3761
3762 return 1;
3763 }
3764
3765 case SOL_ELEMENT:
3766 {
3767 int el1, el2;
3768 mesh->GetTopology().GetSurface2VolumeElement (selnr+1, el1, el2);
3769 el1--;
3770
3771 val = data->data[el1 * data->dist+comp-1];
3772 return 1;
3773 }
3774
3775 case SOL_NONCONTINUOUS:
3776 {
3777 val = 0;
3778 // ?????
3779 return 0;
3780 }
3781
3782 case SOL_SURFACE_ELEMENT:
3783 {
3784 val = data->data[selnr * data->dist + comp-1];
3785 return 1;
3786 }
3787
3788 case SOL_SURFACE_NONCONTINUOUS:
3789 {
3790 const Element2d & el = (*mesh)[selnr];
3791
3792 double lami[8] = { 0.0 };
3793 int np = 0;
3794 val = 0;
3795 int order = data->order;
3796
3797 switch (order)
3798 {
3799 case 0:
3800 return data->data[selnr * data->dist + comp-1];
3801 case 1:
3802 {
3803 switch (el.GetType())
3804 {
3805 case TRIG:
3806 case TRIG6:
3807 {
3808 lami[1] = lam1;
3809 lami[2] = lam2;
3810 lami[0] = 1-lam1-lam2;
3811 np = 3;
3812 break;
3813 }
3814 default:
3815 cerr << "case not impl 234234" << endl;
3816 }
3817 break;
3818 }
3819 case 2:
3820 {
3821 switch (el.GetType())
3822 {
3823 case TRIG:
3824 {
3825 lami[1] = lam1;
3826 lami[2] = lam2;
3827 lami[0] = 1-lam1-lam2;
3828 np = 3;
3829 break;
3830 }
3831 case TRIG6:
3832 {
3833 double lam3 = 1-lam1-lam2;
3834 lami[1] = 2*lam1 * (lam1-0.5);
3835 lami[2] = 2*lam2 * (lam2-0.5);
3836 lami[0] = 2*lam3 * (lam3-0.5);
3837 lami[3] = 4*lam1*lam2;
3838 lami[4] = 4*lam2*lam3;
3839 lami[5] = 4*lam1*lam3;
3840 np = 6;
3841 break;
3842 }
3843 default:
3844 cerr << "case not implented 3234" << endl;
3845 }
3846 break;
3847 }
3848 }
3849
3850 int base;
3851 if (order == 1)
3852 base = 4 * selnr;
3853 else
3854 base = 9 * selnr;
3855
3856 for (int i = 0; i < np; i++)
3857 val += lami[i] * data->data[(base+i) * data->dist + comp-1];
3858
3859 return 1;
3860 }
3861
3862 case SOL_MARKED_ELEMENTS:
3863 {
3864 val = (*mesh)[selnr].TestRefinementFlag();
3865 return 1;
3866 }
3867
3868 case SOL_ELEMENT_ORDER:
3869 {
3870 val = (*mesh)[selnr].GetOrder();
3871 return 1;
3872 }
3873
3874 }
3875 return 0;
3876 }
3877
3878
3879
3880
3881
3882
3883
3884
3885
3886 Vec<3> VisualSceneSolution ::
GetDeformation(ElementIndex elnr,const Point<3> & p) const3887 GetDeformation (ElementIndex elnr, const Point<3> & p) const
3888 {
3889 Vec<3> def;
3890 if (deform && vecfunction != -1)
3891 {
3892 GetValues (soldata[vecfunction], elnr, p(0), p(1), p(2), &def(0));
3893 def *= scaledeform;
3894
3895 if (soldata[vecfunction]->components == 2) def(2) = 0;
3896 }
3897 else
3898 def = 0;
3899 return def;
3900 }
3901
3902
3903 Vec<3> VisualSceneSolution ::
GetSurfDeformation(SurfaceElementIndex elnr,int facetnr,double lam1,double lam2) const3904 GetSurfDeformation (SurfaceElementIndex elnr, int facetnr, double lam1, double lam2) const
3905 {
3906 shared_ptr<Mesh> mesh = GetMesh();
3907
3908 Vec<3> def;
3909 if (deform && vecfunction != -1)
3910 {
3911 // GetSurfValues (soldata[vecfunction], elnr, facetnr, lam1, lam2, &def(0));
3912 double values[6];
3913 GetSurfValues (soldata[vecfunction], elnr, facetnr, lam1, lam2, values);
3914 def = RealVec3d (values, soldata[vecfunction]->iscomplex, imag_part);
3915 def *= scaledeform;
3916
3917 if (soldata[vecfunction]->components == 2) def(2) = 0;
3918 }
3919 else if (deform && scalfunction != -1 && mesh->GetDimension()==2)
3920 { // he: allow for 3d plots of 2d surfaces: usage: turn deformation on
3921 def = 0;
3922 GetSurfValue (soldata[scalfunction], elnr, facetnr, lam1, lam2, scalcomp, def(2));
3923 def *= scaledeform;
3924 }
3925 else
3926 def = 0;
3927 return def;
3928 }
3929
GetPointDeformation(int pnum,Point<3> & p,SurfaceElementIndex elnr) const3930 void VisualSceneSolution :: GetPointDeformation (int pnum, Point<3> & p,
3931 SurfaceElementIndex elnr) const
3932 {
3933 shared_ptr<Mesh> mesh = GetMesh();
3934
3935 p = mesh->Point (pnum+1);
3936 if (deform && vecfunction != -1)
3937 {
3938 const SolData * vsol = soldata[vecfunction];
3939
3940 Vec<3> v(0,0,0);
3941 if (vsol->soltype == SOL_NODAL)
3942 {
3943 v = Vec3d(vsol->data[pnum * vsol->dist],
3944 vsol->data[pnum * vsol->dist+1],
3945 vsol->data[pnum * vsol->dist+2]);
3946 }
3947 else if (vsol->soltype == SOL_SURFACE_NONCONTINUOUS)
3948 {
3949 const Element2d & el = (*mesh)[elnr];
3950 for (int j = 0; j < el.GetNP(); j++)
3951 if (el[j] == pnum+1)
3952 {
3953 int base = (4*elnr+j-1) * vsol->dist;
3954 v = Vec3d(vsol->data[base],
3955 vsol->data[base+1],
3956 vsol->data[base+2]);
3957 }
3958 }
3959
3960 if (vsol->dist == 2) v(2) = 0;
3961
3962 v *= scaledeform;
3963 p += v;
3964 }
3965 }
3966
3967
3968
3969
GetClippingPlaneTrigs(NgArray<ClipPlaneTrig> & trigs,NgArray<ClipPlanePoint> & pts)3970 void VisualSceneSolution :: GetClippingPlaneTrigs (NgArray<ClipPlaneTrig> & trigs,
3971 NgArray<ClipPlanePoint> & pts)
3972 {
3973 shared_ptr<Mesh> mesh = GetMesh();
3974
3975 // static int timer_vals = NgProfiler::CreateTimer ("ClipPlaneTrigs - vertex values");
3976 static int timer1 = NgProfiler::CreateTimer ("ClipPlaneTrigs1");
3977 // static int timer1a = NgProfiler::CreateTimer ("ClipPlaneTrigs1a");
3978 // static int timer2 = NgProfiler::CreateTimer ("ClipPlaneTrigs2");
3979 // static int timer3 = NgProfiler::CreateTimer ("ClipPlaneTrigs3");
3980 // static int timer4 = NgProfiler::CreateTimer ("ClipPlaneTrigs4");
3981 // static int timer4b = NgProfiler::CreateTimer ("ClipPlaneTrigs4b");
3982
3983
3984 NgProfiler::RegionTimer reg1 (timer1);
3985
3986
3987 int ne = mesh->GetNE();
3988
3989 const int edgei[6][2] =
3990 { { 0, 1 }, { 0, 2 }, { 0, 3 },
3991 { 1, 2 }, { 1, 3 }, { 2, 3 } };
3992
3993 double edgelam[6];
3994 // Point<3> edgep[6];
3995 double nodevali[4];
3996
3997 int cntce;
3998 int cpe1 = 0, cpe2 = 0, cpe3 = 0;
3999
4000 // NgArray<Element> loctets;
4001 // NgArray<Element> loctetsloc;
4002 // NgArray<Point<3> > pointsloc;
4003
4004 int n = 1 << subdivisions;
4005 int n3 = (n+1)*(n+1)*(n+1);
4006
4007 NgArray<Point<3> > grid(n3);
4008 NgArray<Point<3> > locgrid(n3);
4009 NgArray<Mat<3,3> > trans(n3);
4010 NgArray<double> val(n3);
4011 NgArray<bool> locposval(n3);
4012 NgArray<int> compress(n3);
4013
4014 // NgProfiler::StartTimer (timer_vals);
4015 NgArray<double,PointIndex::BASE> vertval(mesh->GetNP());
4016 NgArray<bool,PointIndex::BASE> posval(mesh->GetNP());
4017 // for (PointIndex pi = vertval.Begin(); pi < vertval.End(); pi++)
4018 for (PointIndex pi : vertval.Range())
4019 {
4020 Point<3> vert = (*mesh)[pi];
4021 vertval[pi] =
4022 vert(0) * clipplane[0] +
4023 vert(1) * clipplane[1] +
4024 vert(2) * clipplane[2] +
4025 clipplane[3];
4026 posval[pi] = vertval[pi] > 0;
4027 }
4028 // NgProfiler::StopTimer (timer_vals);
4029
4030 INDEX_2_CLOSED_HASHTABLE<int> edges(8*n3); // point nr of edge
4031
4032
4033 for (ElementIndex ei = 0; ei < ne; ei++)
4034 {
4035 // NgProfiler::RegionTimer reg1a (timer1a);
4036 int first_point_of_element = pts.Size();
4037
4038 locgrid.SetSize(n3);
4039 if(vispar.clipdomain > 0 && vispar.clipdomain != (*mesh)[ei].GetIndex()) continue;
4040 if(vispar.donotclipdomain > 0 && vispar.donotclipdomain == (*mesh)[ei].GetIndex()) continue;
4041
4042 ELEMENT_TYPE type = (*mesh)[ei].GetType();
4043 if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID || type == PYRAMID13 || type == PRISM15 || type == HEX20)
4044 {
4045 const Element & el = (*mesh)[ei];
4046
4047 int ii = 0;
4048 int cnt_valid = 0;
4049
4050 // NgProfiler::StartTimer (timer2);
4051
4052
4053 if (!mesh->GetCurvedElements().IsElementHighOrder(ei))
4054 {
4055 bool has_pos = 0, has_neg = 0;
4056
4057 for (int i = 0; i < el.GetNP(); i++)
4058 if (posval[el[i]])
4059 has_pos = 1;
4060 else
4061 has_neg = 1;
4062
4063 if (!has_pos || !has_neg)
4064 {
4065 // NgProfiler::StopTimer (timer2);
4066 continue;
4067 }
4068 }
4069
4070
4071
4072 if (type == TET || type == TET10)
4073 {
4074 for (int ix = 0; ix <= n; ix++)
4075 for (int iy = 0; iy <= n; iy++)
4076 for (int iz = 0; iz <= n; iz++, ii++)
4077 {
4078 if (ix+iy+iz <= n)
4079 {
4080 compress[ii] = cnt_valid;
4081 locgrid[cnt_valid] =
4082 Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
4083 cnt_valid++;
4084 }
4085 else
4086 compress[ii] = -1;
4087 }
4088 }
4089
4090 else
4091
4092 for (int ix = 0; ix <= n; ix++)
4093 for (int iy = 0; iy <= n; iy++)
4094 for (int iz = 0; iz <= n; iz++, ii++)
4095 {
4096 Point<3> ploc;
4097 compress[ii] = ii;
4098
4099 switch (type)
4100 {
4101 case PRISM:
4102 case PRISM12:
4103 case PRISM15:
4104 if (ix+iy <= n)
4105 {
4106 ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
4107 compress[ii] = cnt_valid;
4108 cnt_valid++;
4109 }
4110 else
4111 compress[ii] = -1;
4112 break;
4113 case HEX:
4114 case HEX20:
4115 ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
4116 break;
4117 case PYRAMID:
4118 case PYRAMID13:
4119 ploc = Point<3> (double(ix) / n * (1-double(iz)/n),
4120 double(iy) / n * (1-double(iz)/n),
4121 double(iz)/n);
4122 if (iz == n) ploc = Point<3> (0,0,1-1e-8);
4123 break;
4124 default:
4125 cerr << "clip plane trigs not implemented" << endl;
4126 ploc = Point<3> (0,0,0);
4127 }
4128 if (compress[ii] != -1)
4129 locgrid[compress[ii]] = ploc;
4130 }
4131
4132 if (type != TET && type != TET10 && type != PRISM && type != PRISM12 && type != PRISM15) cnt_valid = n3;
4133
4134 locgrid.SetSize(cnt_valid);
4135
4136 // NgProfiler::StopTimer (timer2);
4137 // NgProfiler::RegionTimer reg4(timer4);
4138
4139 if (mesh->GetCurvedElements().IsHighOrder())
4140 {
4141 // NgProfiler::RegionTimer reg4(timer4);
4142 mesh->GetCurvedElements().
4143 CalcMultiPointElementTransformation (&locgrid, ei, &grid, 0);
4144 }
4145 else
4146 {
4147 // NgProfiler::RegionTimer reg4(timer4b);
4148 Vector shape(el.GetNP());
4149 MatrixFixWidth<3> pointmat(el.GetNP());
4150
4151 for (int k = 0; k < el.GetNP(); k++)
4152 for (int j = 0; j < 3; j++)
4153 pointmat(k,j) = (*mesh)[el[k]](j);
4154
4155 for (int i = 0; i < cnt_valid; i++)
4156 {
4157 el.GetShapeNew<double> (locgrid[i], shape);
4158 Point<3> pglob;
4159 for (int j = 0; j < 3; j++)
4160 {
4161 pglob(j) = 0;
4162 for (int k = 0; k < el.GetNP(); k++)
4163 pglob(j) += shape(k) * pointmat(k,j);
4164 }
4165 grid[i] = pglob;
4166 }
4167 }
4168
4169 // NgProfiler::RegionTimer reg3(timer3);
4170
4171 bool has_pos = false, all_pos = true;
4172
4173 for (int i = 0; i < cnt_valid; i++)
4174 {
4175 val[i] =
4176 grid[i](0) * clipplane[0] +
4177 grid[i](1) * clipplane[1] +
4178 grid[i](2) * clipplane[2] +
4179 clipplane[3];
4180
4181 locposval[i] = val[i] > 0;
4182 has_pos |= locposval[i];
4183 all_pos &= locposval[i];
4184
4185 // if (val[i] > 0) has_pos = 1; else has_neg = 1;
4186 }
4187
4188 // if (!has_pos || !has_neg) continue;
4189 if (!has_pos || all_pos) continue;
4190
4191 edges.DeleteData();
4192
4193 for (int ix = 0; ix < n; ix++)
4194 for (int iy = 0; iy < n; iy++)
4195 for (int iz = 0; iz < n; iz++)
4196 {
4197 int base = iz + (n+1)*iy + (n+1)*(n+1)*ix;
4198 int pi[8] =
4199 { base, base+(n+1)*(n+1), base+(n+1)*(n+1)+(n+1), base+(n+1),
4200 base+1, base+(n+1)*(n+1)+1, base+(n+1)*(n+1)+(n+1)+1, base+(n+1)+1 };
4201
4202 for (int j = 0; j < 8; j++)
4203 pi[j] = compress[pi[j]];
4204
4205 bool has_pos = false, all_pos = true;
4206 for (int j = 0; j < 8; j++)
4207 if (pi[j] != -1)
4208 {
4209 has_pos |= locposval[pi[j]];
4210 all_pos &= locposval[pi[j]];
4211 }
4212 if (!has_pos || all_pos) continue;
4213
4214
4215 const int tets[6][4] =
4216 { { 1, 2, 4, 5 },
4217 { 4, 5, 2, 8 },
4218 { 2, 8, 5, 6 },
4219 { 2, 3, 4, 8 },
4220 { 2, 3, 8, 6 },
4221 { 3, 8, 6, 7 } };
4222
4223 for (int ii = 0; ii < 6; ii++)
4224 {
4225 int teti[4];
4226 for (int k = 0; k < 4; k++)
4227 teti[k] = pi[tets[ii][k]-1];
4228
4229 bool is_valid = true;
4230 for (int j = 0; j < 4; j++)
4231 is_valid &= (teti[j] != -1);
4232 if (!is_valid) continue;
4233
4234 bool has_pos = false, all_pos = true;
4235 for (int j = 0; j < 4; j++)
4236 {
4237 has_pos |= locposval[teti[j]];
4238 all_pos &= locposval[teti[j]];
4239 }
4240 if (!has_pos || all_pos) continue;
4241
4242 for (int j = 0; j < 4; j++)
4243 nodevali[j] = val[teti[j]];
4244
4245 cntce = 0;
4246 for (int j = 0; j < 6; j++)
4247 {
4248 int lpi1 = edgei[j][0];
4249 int lpi2 = edgei[j][1];
4250 if ( (nodevali[lpi1] > 0) != (nodevali[lpi2] > 0) )
4251 {
4252 cntce++;
4253 cpe3 = cpe2;
4254 cpe2 = cpe1;
4255 cpe1 = j;
4256 if (cntce >= 3)
4257 {
4258 ClipPlaneTrig cpt;
4259 cpt.elnr = ei;
4260
4261 for (int k = 0; k < 3; k++)
4262 {
4263 int ednr;
4264 switch (k)
4265 {
4266 case 0: ednr = cpe1; break;
4267 case 1: ednr = cpe2; break;
4268 case 2: ednr = cpe3; break;
4269 }
4270
4271 int pi1 = edgei[ednr][0];
4272 int pi2 = edgei[ednr][1];
4273 int pnr = -1;
4274
4275 INDEX_2 pair (teti[pi1], teti[pi2]);
4276 pair.Sort();
4277 if (edges.Used(pair))
4278 pnr = edges.Get(pair);
4279 else
4280 {
4281 ClipPlanePoint cppt;
4282 cppt.elnr = ei;
4283
4284 edgelam[ednr] = nodevali[pi2] / (nodevali[pi2] - nodevali[pi1]);
4285
4286 Point<3> gp1 = grid[teti[pi1]];
4287 Point<3> gp2 = grid[teti[pi2]];
4288 cppt.p = gp2 + edgelam[ednr] * (gp1-gp2);
4289
4290 Point<3> p1 = locgrid[teti[pi1]];
4291 Point<3> p2 = locgrid[teti[pi2]];
4292 cppt.lami = p2 + edgelam[ednr] * (p1-p2);
4293
4294 pts.Append (cppt);
4295 pnr = pts.Size()-1;
4296 edges.Set (pair, pnr);
4297 }
4298
4299 cpt.points[k].pnr = pnr;
4300 cpt.points[k].locpnr = pnr-first_point_of_element;
4301 }
4302
4303 trigs.Append (cpt);
4304 }
4305 }
4306 }
4307 }
4308 }
4309 }
4310
4311 else
4312 { // other elements not supported (JS, June 2007)
4313 continue; // return;
4314 }
4315
4316 }
4317 }
4318
GetClippingPlaneGrid(NgArray<ClipPlanePoint> & pts)4319 void VisualSceneSolution :: GetClippingPlaneGrid (NgArray<ClipPlanePoint> & pts)
4320 {
4321 shared_ptr<Mesh> mesh = GetMesh();
4322
4323 Vec3d n(clipplane[0], clipplane[1], clipplane[2]);
4324
4325 double mu = -clipplane[3] / n.Length2();
4326 Point3d p(mu*n.X(), mu * n.Y(), mu * n.Z());
4327
4328 // n /= n.Length();
4329 n.Normalize();
4330 Vec3d t1, t2;
4331 n.GetNormal (t1);
4332 t2 = Cross (n, t1);
4333
4334 double xi1, xi2;
4335
4336 double xi1mid = (center - p) * t1;
4337 double xi2mid = (center - p) * t2;
4338
4339 pts.SetSize(0);
4340
4341 for (xi1 = xi1mid-rad+xoffset/gridsize; xi1 <= xi1mid+rad+xoffset/gridsize; xi1 += rad / gridsize)
4342 for (xi2 = xi2mid-rad+yoffset/gridsize; xi2 <= xi2mid+rad+yoffset/gridsize; xi2 += rad / gridsize)
4343 {
4344 Point3d hp = p + xi1 * t1 + xi2 * t2;
4345
4346 int cindex(-1);
4347 bool allowindex(true);
4348 if(vispar.clipdomain > 0)
4349 {
4350 cindex = vispar.clipdomain;
4351 }
4352 else if(vispar.donotclipdomain > 0)
4353 {
4354 allowindex = false;
4355 cindex = vispar.donotclipdomain;
4356 }
4357
4358 double lami[3];
4359 int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex)-1;
4360
4361 if (elnr != -1)
4362 {
4363 ClipPlanePoint cpp;
4364 cpp.p = hp;
4365 cpp.elnr = elnr;
4366 cpp.lami(0) = lami[0];
4367 cpp.lami(1) = lami[1];
4368 cpp.lami(2) = lami[2];
4369 pts.Append (cpp);
4370 }
4371 }
4372 };
4373
4374
4375
4376
DrawClipPlaneTrigs()4377 void VisualSceneSolution :: DrawClipPlaneTrigs ()
4378 {
4379 shared_ptr<Mesh> mesh = GetMesh();
4380
4381 #ifdef PARALLELGL
4382
4383 if (id == 0 && ntasks > 1)
4384 {
4385 InitParallelGL();
4386
4387 NgArray<int> parlists (ntasks);
4388
4389 MyMPI_SendCmd ("redraw");
4390 MyMPI_SendCmd ("clipplanetrigs");
4391
4392 for ( int dest = 1; dest < ntasks; dest++ )
4393 MyMPI_Recv (parlists[dest], dest, MPI_TAG_VIS);
4394
4395 if (clipplanelist_scal)
4396 glDeleteLists (clipplanelist_scal, 1);
4397
4398 clipplanelist_scal = glGenLists (1);
4399 glNewList (clipplanelist_scal, GL_COMPILE);
4400
4401 for ( int dest = 1; dest < ntasks; dest++ )
4402 glCallList (parlists[dest]);
4403
4404 glEndList();
4405 return;
4406 }
4407 #endif
4408
4409
4410
4411
4412
4413 if (clipplanelist_scal)
4414 glDeleteLists (clipplanelist_scal, 1);
4415
4416 clipplanelist_scal = glGenLists (1);
4417 glNewList (clipplanelist_scal, GL_COMPILE);
4418
4419
4420 NgArray<ClipPlaneTrig> trigs;
4421 NgArray<ClipPlanePoint> points;
4422 GetClippingPlaneTrigs (trigs, points);
4423
4424 glNormal3d (-clipplane[0], -clipplane[1], -clipplane[2]);
4425 glColor3d (1.0, 1.0, 1.0);
4426
4427 SetTextureMode (usetexture);
4428
4429 SolData * sol = NULL;
4430
4431 if (scalfunction != -1)
4432 sol = soldata[scalfunction];
4433
4434 if (sol -> draw_volume)
4435 {
4436 glBegin (GL_TRIANGLES);
4437
4438 int maxlpnr = 0;
4439 for (int i = 0; i < trigs.Size(); i++)
4440 for (int j = 0; j < 3; j++)
4441 maxlpnr = max2 (maxlpnr, trigs[i].points[j].locpnr);
4442
4443 NgArray<double> vals(maxlpnr+1);
4444 NgArray<complex<double> > valsc(maxlpnr+1);
4445 NgArray<int> elnrs(maxlpnr+1);
4446 NgArray<bool> trigok(maxlpnr+1);
4447 NgArray<Point<3> > locpoints(maxlpnr+1);
4448 NgArray<Point<3> > globpoints(maxlpnr+1);
4449 NgArray<Mat<3> > jacobi(maxlpnr+1);
4450 NgArray<double> mvalues( (maxlpnr+1) * sol->components);
4451 trigok = false;
4452 elnrs = -1;
4453
4454 Point<3> p[3];
4455 // double val[3];
4456 // complex<double> valc[3];
4457 int lastelnr = -1;
4458 int nlp = -1;
4459 bool ok = false;
4460
4461 for (int i = 0; i < trigs.Size(); i++)
4462 {
4463 const ClipPlaneTrig & trig = trigs[i];
4464 if (trig.elnr != lastelnr)
4465 {
4466 lastelnr = trig.elnr;
4467 nlp = -1;
4468
4469 for (int ii = i; ii < trigs.Size(); ii++)
4470 {
4471 if (trigs[ii].elnr != trig.elnr) break;
4472 for (int j = 0; j < 3; j++)
4473 nlp = max (nlp, trigs[ii].points[j].locpnr);
4474 }
4475 nlp++;
4476 locpoints.SetSize (nlp);
4477
4478 for (int ii = i; ii < trigs.Size(); ii++)
4479 {
4480 if (trigs[ii].elnr != trig.elnr) break;
4481 for (int j = 0; j < 3; j++)
4482 locpoints[trigs[ii].points[j].locpnr] = points[trigs[ii].points[j].pnr].lami;
4483 }
4484
4485 mesh->GetCurvedElements().
4486 CalcMultiPointElementTransformation (&locpoints, trig.elnr,
4487 &globpoints, &jacobi);
4488
4489 bool
4490 drawelem = GetMultiValues (sol, trig.elnr, -1, nlp,
4491 &locpoints[0](0), &locpoints[1](0)-&locpoints[0](0),
4492 &globpoints[0](0), &globpoints[1](0)-&globpoints[0](0),
4493 &jacobi[0](0), &jacobi[1](0)-&jacobi[0](0),
4494 &mvalues[0], sol->components);
4495
4496 // cout << "have multivalues, comps = " << sol->components << endl;
4497
4498 // if (!drawelem) ok = false;
4499 ok = drawelem;
4500 if (usetexture != 2 || !sol->iscomplex)
4501 for (int ii = 0; ii < nlp; ii++)
4502 vals[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]);
4503 else
4504 for (int ii = 0; ii < nlp; ii++)
4505 valsc[ii] = complex<double> (mvalues[ii*sol->components + scalcomp-1],
4506 mvalues[ii*sol->components + scalcomp]);
4507 }
4508
4509 if(ok)
4510 for(int j=0; j<3; j++)
4511 {
4512 if (usetexture != 2 || !sol->iscomplex)
4513 SetOpenGlColor (vals[trig.points[j].locpnr]);
4514 else
4515 glTexCoord2f ( valsc[trig.points[j].locpnr].real(),
4516 valsc[trig.points[j].locpnr].imag() );
4517
4518 p[j] = points[trig.points[j].pnr].p;
4519
4520 if (deform)
4521 {
4522 Point<3> ploc = points[trig.points[j].pnr].lami;
4523 p[j] += GetDeformation (trig.elnr, ploc);
4524 }
4525
4526 glVertex3dv (p[j]);
4527 }
4528
4529 }
4530 glEnd();
4531 }
4532 glEndList ();
4533
4534
4535 #ifdef PARALLELGL
4536 glFinish();
4537 if (id > 0)
4538 MyMPI_Send (clipplanelist_scal, 0, MPI_TAG_VIS);
4539 #endif
4540 }
4541
4542
4543
4544
4545
4546
4547
4548
4549
4550
4551 void VisualSceneSolution ::
SetOpenGlColor(double val)4552 SetOpenGlColor(double val)
4553 {
4554 if (usetexture == 1 && !logscale)
4555 {
4556 glTexCoord1f ( val );
4557 return;
4558 }
4559
4560 double valmin = minval;
4561 double valmax = maxval;
4562
4563 double value;
4564
4565 if (!logscale)
4566 value = (val - valmin) / (valmax - valmin);
4567 else
4568 {
4569 if (valmax <= 0) valmax = 1;
4570 if (valmin <= 0) valmin = 1e-4 * valmax;
4571 value = (log(fabs(val)) - log(valmin)) / (log(valmax) - log(valmin));
4572 }
4573
4574 if (!invcolor)
4575 value = 1 - value;
4576
4577
4578 if (value > 1) value = 1;
4579 if (value < 0) value = 0;
4580
4581 value *= 4;
4582
4583 static const double colp[][3] =
4584 {
4585 { 1, 0, 0 },
4586 { 1, 1, 0 },
4587 { 0, 1, 0 },
4588 { 0, 1, 1 },
4589 { 0, 0, 1 },
4590 { 1, 0, 1 },
4591 { 1, 0, 0 },
4592 };
4593
4594 int i = int(value);
4595 double r = value - i;
4596
4597 GLdouble col[3];
4598 for (int j = 0; j < 3; j++)
4599 col[j] = (1-r) * colp[i][j] + r * colp[i+1][j];
4600
4601 glColor3dv (col);
4602 }
4603
4604
4605
4606 void VisualSceneSolution ::
SetTextureMode(int texturemode) const4607 SetTextureMode (int texturemode) const
4608 {
4609 switch (texturemode)
4610 {
4611 case 0:
4612 glDisable (GL_TEXTURE_1D);
4613 glDisable (GL_TEXTURE_2D);
4614 break;
4615 case 1:
4616 glEnable (GL_TEXTURE_1D);
4617 glDisable (GL_TEXTURE_2D);
4618 glColor3d (1.0, 1.0, 1.0);
4619 break;
4620 case 2:
4621 glDisable (GL_TEXTURE_1D);
4622 glEnable (GL_TEXTURE_2D);
4623 glColor3d (1.0, 1.0, 1.0);
4624 break;
4625 }
4626 }
4627
4628
4629
4630
4631 void VisualSceneSolution ::
DrawCone(const Point<3> & p1,const Point<3> & p2,double r)4632 DrawCone (const Point<3> & p1, const Point<3> & p2, double r)
4633 {
4634 int n = 10, i;
4635 Vec<3> p1p2 = p2 - p1;
4636
4637 p1p2.Normalize();
4638 Vec<3> p2p1 = -p1p2;
4639
4640 Vec<3> t1 = p1p2.GetNormal();
4641 Vec<3> t2 = Cross (p1p2, t1);
4642
4643 Point<3> oldp = p1 + r * t1;
4644 Vec<3> oldn = t1;
4645
4646 Point<3> p;
4647 Vec<3> normal;
4648
4649 Mat<2> rotmat;
4650 Vec<2> cs, newcs;
4651 cs(0) = 1;
4652 cs(1) = 0;
4653 rotmat(0,0) = rotmat(1,1) = cos(2*M_PI/n);
4654 rotmat(1,0) = sin(2*M_PI/n);
4655 rotmat(0,1) = -rotmat(1,0);
4656
4657 glBegin (GL_TRIANGLES);
4658 for (i = 1; i <= n; i++)
4659 {
4660 /*
4661 phi = 2 * M_PI * i / n;
4662 normal = cos(phi) * t1 + sin(phi) * t2;
4663 */
4664 newcs = rotmat * cs;
4665 cs = newcs;
4666 normal = cs(0) * t1 + cs(1) * t2;
4667
4668 p = p1 + r * normal;
4669
4670 // cone
4671 glNormal3dv (normal);
4672 glVertex3dv (p);
4673 glVertex3dv (p2);
4674 glNormal3dv (oldn);
4675 glVertex3dv (oldp);
4676
4677 // base-circle
4678 glNormal3dv (p2p1);
4679 glVertex3dv (p);
4680 glVertex3dv (p1);
4681 glVertex3dv (oldp);
4682
4683 oldp = p;
4684 oldn = normal;
4685 }
4686 glEnd ();
4687 }
4688
4689
4690
4691 void VisualSceneSolution ::
DrawCylinder(const Point<3> & p1,const Point<3> & p2,double r)4692 DrawCylinder (const Point<3> & p1, const Point<3> & p2, double r)
4693 {
4694 int n = 10, i;
4695 Vec<3> p1p2 = p2 - p1;
4696
4697 p1p2.Normalize();
4698 // Vec<3> p2p1 = -p1p2;
4699
4700 Vec<3> t1 = p1p2.GetNormal();
4701 Vec<3> t2 = Cross (p1p2, t1);
4702
4703 Point<3> oldhp1 = p1 + r * t1;
4704 Point<3> oldhp2 = p2 + r * t1;
4705 Vec<3> oldn = t1;
4706
4707 Point<3> hp1, hp2;
4708 Vec<3> normal;
4709
4710 Mat<2> rotmat;
4711 Vec<2> cs, newcs;
4712 cs(0) = 1;
4713 cs(1) = 0;
4714 rotmat(0,0) = rotmat(1,1) = cos(2*M_PI/n);
4715 rotmat(1,0) = sin(2*M_PI/n);
4716 rotmat(0,1) = -rotmat(1,0);
4717
4718 glBegin (GL_QUADS);
4719 for (i = 1; i <= n; i++)
4720 {
4721 newcs = rotmat * cs;
4722 cs = newcs;
4723 normal = cs(0) * t1 + cs(1) * t2;
4724
4725 hp1 = p1 + r * normal;
4726 hp2 = p2 + r * normal;
4727
4728 // cylinder
4729 glNormal3dv (normal);
4730
4731 glVertex3dv (hp1);
4732 glVertex3dv (hp2);
4733 glVertex3dv (oldhp2);
4734 glVertex3dv (oldhp1);
4735
4736 oldhp1 = hp1;
4737 oldhp2 = hp2;
4738 oldn = normal;
4739 }
4740 glEnd ();
4741 }
4742
4743
4744
4745
4746
4747
4748
4749
4750
4751
4752
4753
4754
MouseDblClick(int px,int py)4755 void VisualSceneSolution :: MouseDblClick (int px, int py)
4756 {
4757 auto mesh = GetMesh();
4758 auto dim = mesh->GetDimension();
4759
4760 auto formatComplex = [](double real, double imag)
4761 {
4762 return ToString(real) + (imag < 0 ? "" : "+") + ToString(imag) + "j";
4763 };
4764
4765 auto printScalValue = [&formatComplex]
4766 (SolData & sol, int comp, double value, double imag=0., bool iscomplex=false)
4767 {
4768 if(sol.components>1)
4769 {
4770 if(comp==0)
4771 cout << "func(" << sol.name << ")";
4772 else
4773 cout << sol.name << "["+ToString(comp)+"]";
4774 }
4775 else
4776 cout << sol.name;
4777 cout << " = " << (iscomplex ? formatComplex(value, imag) : ToString(value)) << endl;
4778 };
4779
4780 auto printVecValue = [&formatComplex]
4781 (SolData & sol, FlatArray<double> values)
4782 {
4783 if(sol.iscomplex)
4784 {
4785 cout << sol.name << " = ( " << formatComplex(values[0], values[1]);
4786 for(int i = 2; i < values.Size(); i+=2)
4787 cout << ", " << formatComplex(values[i], values[i+1]);
4788 cout << " )" << endl;
4789 }
4790 else
4791 {
4792 cout << sol.name << " = ( " << values[0];
4793 for(int i = 1; i < values.Size(); i++)
4794 cout << ", " << values[i];
4795 cout << " )" << endl;
4796 }
4797 };
4798
4799 // Check if clipping plane is drawn at current mouse cursor position
4800 if(dim==3 && clipsolution && vispar.clipping.enable)
4801 {
4802 GLint viewport[4];
4803 GLdouble projection[16];
4804 glGetDoublev(GL_PROJECTION_MATRIX, &projection[0]);
4805
4806 glGetIntegerv(GL_VIEWPORT, &viewport[0]);
4807
4808 int hy = viewport[3]-py;
4809
4810 // manually intersect the view vector with the clipping plane (also working if clipping vectors are shown)
4811 Point<3> p_clipping_plane;
4812 gluUnProject(px, hy, 1.0, transformationmat, projection, viewport,
4813 &p_clipping_plane[0], &p_clipping_plane[1], &p_clipping_plane[2]);
4814
4815 Point<3> eye;
4816 gluUnProject( (viewport[2]-viewport[0])/2 , (viewport[3]-viewport[1])/2,
4817 0.0, transformationmat, projection, viewport, &eye[0], &eye[1], &eye[2]);
4818
4819 Vec<3> n{vispar.clipping.normal};
4820 n.Normalize();
4821 Vec<3> view = p_clipping_plane-eye;
4822
4823 // check if we look at the clipping plane from the right direction
4824 if(n*view > 1e-8)
4825 {
4826 double lam = vispar.clipping.dist - Vec<3>{eye}*n;
4827 lam /= n*view;
4828 p_clipping_plane = eye + lam*view;
4829
4830 double lami[3];
4831 if(auto el3d = mesh->GetElementOfPoint( p_clipping_plane, lami ))
4832 {
4833 cout << endl << "Selected point " << p_clipping_plane << " on clipping plane" << endl;
4834
4835 bool have_scal_func = scalfunction!=-1 && soldata[scalfunction]->draw_volume;
4836 bool have_vec_func = vecfunction!=-1 && soldata[vecfunction]->draw_volume;
4837
4838 if(have_scal_func)
4839 {
4840 auto & sol = *soldata[scalfunction];
4841 double val;
4842 double imag = 0;
4843 int rcomponent = scalcomp;
4844 int comp = scalcomp;
4845 if(sol.iscomplex && rcomponent != 0)
4846 {
4847 rcomponent = 2 * ((rcomponent-1)/2) + 1;
4848 GetValue(&sol, el3d-1, lami[0], lami[1], lami[2], rcomponent+1,
4849 imag);
4850 comp = (scalcomp-1)/2 + 1;
4851 }
4852 GetValue(&sol, el3d-1, lami[0], lami[1], lami[2], rcomponent, val);
4853 printScalValue(sol, comp, val, imag, sol.iscomplex && comp > 0);
4854 }
4855 if(vecfunction!=-1 && soldata[vecfunction]->draw_volume)
4856 {
4857 auto & sol = *soldata[vecfunction];
4858 ArrayMem<double, 10> values(sol.components);
4859 GetValues(&sol, el3d-1, lami[0], lami[1], lami[2], &values[0]);
4860 printVecValue(sol, values);
4861 }
4862 return;
4863 }
4864 }
4865 }
4866
4867 // no point on clipping plane found -> continue searching for surface element
4868
4869 Point<3> p;
4870 bool found_point = vsmesh.Unproject(px, py, p);
4871 if(!found_point)
4872 return;
4873
4874 if(selelement==0)
4875 return;
4876
4877 double lami[3] = {0.0, 0.0, 0.0};
4878 // Check if unprojected Point is close to surface element (eps of 1e-3 due to z-Buffer accuracy)
4879 bool found_2del = false;
4880 if(mesh->PointContainedIn2DElement(p, lami, selelement, false && fabs(lami[2])<1e-3))
4881 {
4882 // Found it, use coordinates of point projected to surface element
4883 mesh->GetCurvedElements().CalcSurfaceTransformation({1.0-lami[0]-lami[1], lami[0]}, selelement-1, p);
4884 found_2del = true;
4885 }
4886 cout << endl << "Selected point " << p << " on surface" << endl;
4887
4888 if(!found_2del)
4889 return;
4890
4891 bool have_scal_func = scalfunction!=-1 && soldata[scalfunction]->draw_surface;
4892 bool have_vec_func = vecfunction!=-1 && soldata[vecfunction]->draw_surface;
4893
4894 if(have_scal_func)
4895 {
4896 auto & sol = *soldata[scalfunction];
4897 double val;
4898 double imag = 0;
4899 int rcomponent = scalcomp;
4900 int comp = scalcomp;
4901 if(sol.iscomplex && rcomponent != 0)
4902 {
4903 rcomponent = 2 * ((rcomponent-1)/2) + 1;
4904 GetSurfValue(&sol, selelement-1, -1, 1.0-lami[0]-lami[1], lami[0], rcomponent+1, imag);
4905 comp = (scalcomp-1)/2 + 1;
4906 }
4907 GetSurfValue(&sol, selelement-1, -1, 1.0-lami[0]-lami[1], lami[0], rcomponent, val);
4908 printScalValue(sol, comp, val, imag, sol.iscomplex && comp > 0);
4909 }
4910 if(have_vec_func)
4911 {
4912 auto & sol = *soldata[vecfunction];
4913 ArrayMem<double, 10> values(sol.components);
4914 GetSurfValues(&sol, selelement-1, -1, 1.0-lami[0]-lami[1], lami[0], &values[0]);
4915 printVecValue(sol, values);
4916 }
4917 }
4918
4919
4920
4921 #ifdef PARALLELGL
4922
Broadcast()4923 void VisualSceneSolution :: Broadcast ()
4924 {
4925 MPI_Datatype type;
4926 int blocklen[] =
4927 {
4928 1, 1, 1, 1,
4929 1, 1, 1, 1,
4930 1, 1, 1, 1,
4931 1, 4, 1, 1,
4932 1
4933 };
4934 MPI_Aint displ[] = { (char*)&usetexture - (char*)this,
4935 (char*)&clipsolution - (char*)this,
4936 (char*)&scalfunction - (char*)this,
4937 (char*)&scalcomp - (char*)this,
4938
4939 (char*)&vecfunction - (char*)this,
4940 (char*)&gridsize - (char*)this,
4941 (char*)&autoscale - (char*)this,
4942 (char*)&logscale - (char*)this,
4943
4944 (char*)&minval - (char*)this,
4945 (char*)&maxval - (char*)this,
4946 (char*)&numisolines - (char*)this,
4947 (char*)&subdivisions - (char*)this,
4948
4949 (char*)&evalfunc - (char*)this,
4950 (char*)&clipplane[0] - (char*)this,
4951 (char*)&multidimcomponent - (char*)this,
4952 (char*)&deform - (char*)this,
4953
4954 (char*)&scaledeform - (char*)this
4955 };
4956
4957
4958 MPI_Datatype types[] = {
4959 MPI_INT, MPI_INT, MPI_INT, MPI_INT,
4960 MPI_INT, MPI_INT, MPI_INT, MPI_INT,
4961 MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT,
4962 MPI_INT, MPI_DOUBLE, MPI_INT, MPI_INT,
4963 MPI_DOUBLE
4964 };
4965
4966 MPI_Type_create_struct (17, blocklen, displ, types, &type);
4967 MPI_Type_commit ( &type );
4968
4969 MPI_Bcast (this, 1, type, 0, MPI_COMM_WORLD);
4970 MPI_Type_free (&type);
4971 }
4972
4973 #endif
4974
4975 }
4976
4977
4978 #include "../include/nginterface.h"
4979
Ng_ClearSolutionData()4980 void Ng_ClearSolutionData ()
4981 {
4982 #ifdef OPENGL
4983 // if (nodisplay) return;
4984 // netgen::vssolution.ClearSolutionData();
4985 netgen::GetVSSolution().ClearSolutionData();
4986 #endif
4987 }
4988
Ng_InitSolutionData(Ng_SolutionData * soldata)4989 void Ng_InitSolutionData (Ng_SolutionData * soldata)
4990 {
4991 // soldata -> name = NULL;
4992 soldata -> data = NULL;
4993 soldata -> components = 1;
4994 soldata -> dist = 1;
4995 soldata -> order = 1;
4996 soldata -> iscomplex = 0;
4997 soldata -> draw_surface = 1;
4998 soldata -> draw_volume = 1;
4999 soldata -> soltype = NG_SOLUTION_NODAL;
5000 soldata -> solclass = 0;
5001 }
5002
Ng_SetSolutionData(Ng_SolutionData * soldata)5003 void Ng_SetSolutionData (Ng_SolutionData * soldata)
5004 {
5005 #ifdef OPENGL
5006 // if (nodisplay) return;
5007 // vssolution.ClearSolutionData ();
5008 netgen::VisualSceneSolution::SolData * vss = new netgen::VisualSceneSolution::SolData;
5009
5010 // vss->name = new char[strlen (soldata->name)+1];
5011 // strcpy (vss->name, soldata->name);
5012 vss->name = soldata->name;
5013 vss->data = soldata->data;
5014 vss->components = soldata->components;
5015 vss->dist = soldata->dist;
5016 vss->order = soldata->order;
5017 vss->iscomplex = bool(soldata->iscomplex);
5018 vss->draw_surface = soldata->draw_surface;
5019 vss->draw_volume = soldata->draw_volume;
5020 vss->soltype = netgen::VisualSceneSolution::SolType (soldata->soltype);
5021 vss->solclass = soldata->solclass;
5022 // netgen::vssolution.AddSolutionData (vss);
5023 netgen::GetVSSolution().AddSolutionData (vss);
5024 #endif
5025 }
5026
5027
5028
5029 namespace netgen
5030 {
5031 extern void Render (bool blocking);
5032 }
5033
Ng_Redraw(bool blocking)5034 void Ng_Redraw (bool blocking)
5035 {
5036 #ifdef OPENGL
5037 //netgen::vssolution.UpdateSolutionTimeStamp();
5038 netgen::GetVSSolution().UpdateSolutionTimeStamp();
5039 netgen::Render(blocking);
5040 #endif
5041 }
5042
5043 #ifdef OPENGL
5044 #ifdef WIN32
5045 void (*glBindBuffer) (GLenum a, GLuint b);
5046 void (*glDeleteBuffers) (GLsizei a, const GLuint *b);
5047 void (*glGenBuffers) (GLsizei a, GLuint *b);
5048 void (*glBufferData) (GLenum a, GLsizeiptr b, const GLvoid *c, GLenum d);
5049 void (*glBufferSubData) (GLenum a, GLintptr b, GLsizeiptr c, const GLvoid *d);
5050
5051 GLenum (*glCheckFramebufferStatus) (GLenum target);
5052 void (*glBindFramebuffer) (GLenum target, GLuint framebuffer);
5053 void (*glBindRenderbuffer) (GLenum target, GLuint renderbuffer);
5054 void (*glDeleteFramebuffers) (GLsizei n, const GLuint *framebuffers);
5055 void (*glDeleteRenderbuffers) (GLsizei n, const GLuint *renderbuffers);
5056 void (*glGenFramebuffers) (GLsizei n, GLuint *framebuffers);
5057 void (*glGenRenderbuffers) (GLsizei n, GLuint *renderbuffers);
5058 void (*glRenderbufferStorage) (GLenum target, GLenum internalformat, GLsizei width, GLsizei height);
5059 void (*glFramebufferRenderbuffer) (GLenum target, GLenum attachment, GLenum renderbuffertarget, GLuint renderbuffer);
5060
LoadOpenGLFunctionPointers()5061 DLL_HEADER void LoadOpenGLFunctionPointers() {
5062 #ifdef USE_BUFFERS
5063 glBindBuffer = (decltype(glBindBuffer)) wglGetProcAddress("glBindBuffer");
5064 glBufferSubData = (decltype(glBufferSubData)) wglGetProcAddress("glBufferSubData");
5065 glBufferData = (decltype(glBufferData)) wglGetProcAddress("glBufferData");
5066 glDeleteBuffers = (decltype(glDeleteBuffers)) wglGetProcAddress("glDeleteBuffers");
5067 glGenBuffers = (decltype(glGenBuffers)) wglGetProcAddress("glGenBuffers");
5068 if(!glBindBuffer) throw std::runtime_error("Could not load OpenGL functions!");
5069 #endif
5070
5071 glCheckFramebufferStatus = (decltype(glCheckFramebufferStatus )) wglGetProcAddress("glCheckFramebufferStatus");
5072 glBindFramebuffer = (decltype(glBindFramebuffer )) wglGetProcAddress("glBindFramebuffer");
5073 glBindRenderbuffer = (decltype(glBindRenderbuffer )) wglGetProcAddress("glBindRenderbuffer");
5074 glDeleteFramebuffers = (decltype(glDeleteFramebuffers )) wglGetProcAddress("glDeleteFramebuffers");
5075 glDeleteRenderbuffers = (decltype(glDeleteRenderbuffers )) wglGetProcAddress("glDeleteRenderbuffers");
5076 glGenFramebuffers = (decltype(glGenFramebuffers )) wglGetProcAddress("glGenFramebuffers");
5077 glGenRenderbuffers = (decltype(glGenRenderbuffers )) wglGetProcAddress("glGenRenderbuffers");
5078 glRenderbufferStorage = (decltype(glRenderbufferStorage )) wglGetProcAddress("glRenderbufferStorage");
5079 glFramebufferRenderbuffer = (decltype(glFramebufferRenderbuffer )) wglGetProcAddress("glFramebufferRenderbuffer");
5080 }
5081 #else // WIN32
LoadOpenGLFunctionPointers()5082 DLL_HEADER void LoadOpenGLFunctionPointers() { }
5083 #endif // WIN32
5084 #endif // OPENGL
5085