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