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