src/PIPELINE/VISU_Streamer.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "VISU_Streamer.hxx"
  21
  22 #include "vtkCell.h"
  23 #include "vtkDataSet.h"
  24 #include "vtkDoubleArray.h"
  25 #include "vtkExecutive.h"
  26 #include "vtkGenericCell.h"
  27 #include "vtkInformation.h"
  28 #include "vtkInterpolatedVelocityField.h"
  29 #include "vtkMath.h"
  30 #include "vtkMultiThreader.h"
  31 #include "vtkObjectFactory.h"
  32 #include "vtkPointData.h"
  33 #include "vtkRungeKutta2.h"
  34
  35 vtkCxxRevisionMacro(VISU_Streamer, "$Revision$");
  36
  37 #define VTK_START_FROM_POSITION 0
  38 #define VTK_START_FROM_LOCATION 1
  39
  40 struct VISU_StreamerThreadStruct
  41 {
  42   VISU_Streamer *Filter;
  43   vtkDataSet *Input;
  44   vtkDataSet *Source;
  45 };
  46
  47 vtkStreamer::StreamArray::StreamArray()
  48 {
  49   this->MaxId = -1;
  50   this->Array = new vtkStreamer::StreamPoint[1000];
  51   this->Size = 1000;
  52   this->Extend = 5000;
  53   this->Direction = VTK_INTEGRATE_FORWARD;
  54 }
  55
  56 vtkStreamer::StreamPoint *vtkStreamer::StreamArray::Resize(vtkIdType sz)
  57 {
  58   vtkStreamer::StreamPoint *newArray;
  59   vtkIdType newSize;
  60
  61   if (sz >= this->Size)
  62     {
  63     newSize = this->Size +
  64       this->Extend*(((sz-this->Size)/this->Extend)+1);
  65     }
  66   else
  67     {
  68     newSize = sz;
  69     }
  70
  71   newArray = new vtkStreamer::StreamPoint[newSize];
  72
  73   memcpy(newArray, this->Array,
  74          (sz < this->Size ? sz : this->Size) * sizeof(vtkStreamer::StreamPoint));
  75
  76   this->Size = newSize;
  77   delete [] this->Array;
  78   this->Array = newArray;
  79
  80   return this->Array;
  81 }
  82
  83 // Construct object to start from position (0,0,0); integrate forward; terminal
  84 // speed 0.0; vorticity computation off; integrations step length 0.2; and
  85 // maximum propagation time 100.0.
  86 VISU_Streamer::VISU_Streamer()
  87 {
  88 }
  89
  90 VISU_Streamer::~VISU_Streamer()
  91 {
  92 }
  93
  94 static const double VTK_EPSILON=1E-12;
  95
  96 VTK_THREAD_RETURN_TYPE VISU_Streamer::ThreadedIntegrate( void *arg )
  97 {
  98   VISU_Streamer              *self;
  99   VISU_StreamerThreadStruct  *str;
 100   int                      thread_count;
 101   int                      thread_id;
 102   vtkStreamer::StreamArray *streamer;
 103   vtkStreamer::StreamPoint *sNext = 0, *sPtr;
 104   vtkStreamer::StreamPoint pt1, pt2;
 105   int                      i;
 106   vtkIdType                idxNext, ptId;
 107   double                   d, step, dir;
 108   double                   xNext[3], vel[3];
 109   double                   *cellVel;
 110   double                   derivs[9];
 111   double                   *w, pcoords[3];
 112   double                   coords[4];
 113   vtkDataSet               *input;
 114   vtkGenericCell           *cell;
 115   vtkPointData             *pd;
 116   vtkDataArray             *inScalars;
 117   vtkDataArray             *inVectors;
 118   vtkDoubleArray           *cellVectors;
 119   vtkDataArray             *cellScalars=0;
 120   double tOffset, vort[3];
 121   double err;
 122   int counter=0;
 123
 124   thread_id = ((vtkMultiThreader::ThreadInfo *)(arg))->ThreadID;
 125   thread_count = ((vtkMultiThreader::ThreadInfo *)(arg))->NumberOfThreads;
 126   str = (VISU_StreamerThreadStruct*)(((vtkMultiThreader::ThreadInfo *)(arg))->UserData);
 127   self = str->Filter;
 128
 129   input     = str->Input;
 130   pd        = input->GetPointData();
 131   inScalars = pd->GetScalars();
 132   inVectors = pd->GetVectors();
 133
 134   cell = vtkGenericCell::New();
 135   cellVectors = vtkDoubleArray::New();
 136   cellVectors->SetNumberOfComponents(3);
 137   cellVectors->Allocate(3*VTK_CELL_SIZE);
 138   if (inScalars)
 139     {
 140     cellScalars = inScalars->NewInstance();
 141     cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
 142     cellScalars->Allocate(inScalars->GetNumberOfComponents()*VTK_CELL_SIZE);
 143     }
 144
 145   w = new double[input->GetMaxCellSize()];
 146
 147   // Set the function set to be integrated
 148   vtkInterpolatedVelocityField* func = vtkInterpolatedVelocityField::New();
 149   func->AddDataSet(input);
 150
 151   if (self->GetIntegrator() == 0)
 152     {
 153     vtkGenericWarningMacro("No integrator is specified.");
 154     return VTK_THREAD_RETURN_VALUE;
 155     }
 156
 157   // Create a new integrator, the type is the same as Integrator
 158   vtkInitialValueProblemSolver* integrator =
 159     self->GetIntegrator()->NewInstance();
 160   integrator->SetFunctionSet(func);
 161
 162   // Used to avoid calling these function many times during
 163   // the integration
 164   double termspeed = self->GetTerminalSpeed();
 165   double maxtime = self->GetMaximumPropagationTime();
 166   double savePointInterval = self->GetSavePointInterval();
 167
 168   // For each streamer, integrate in appropriate direction
 169   // Do only the streamers that this thread should handle.
 170   for (ptId=0; ptId < self->GetNumberOfStreamers(); ptId++)
 171     {
 172     if ( ptId % thread_count == thread_id )
 173       {
 174       // Get starting step
 175       streamer = self->GetStreamers() + ptId;
 176       sPtr = streamer->GetStreamPoint(0);
 177       if ( sPtr->cellId < 0 )
 178         {
 179         continue;
 180         }
 181       // Set the last cell id in the vtkInterpolatedVelocityField
 182       // object to speed up FindCell calls
 183       func->SetLastCellId(sPtr->cellId);
 184
 185       dir = streamer->Direction;
 186
 187       // Copy the first point
 188       pt1 = *sPtr;
 189       pt2 = *sPtr;
 190       tOffset = pt1.t;
 191
 192       //integrate until time has been exceeded
 193       while ( pt1.cellId >= 0 && pt1.speed > termspeed && pt1.t <  maxtime )
 194         {
 195
 196         if ( counter++ % 1000 == 0 )
 197           {
 198           if (!thread_id)
 199             {
 200             self->UpdateProgress((double)ptId/self->GetNumberOfStreamers()
 201                                  +pt1.t/maxtime/self->GetNumberOfStreamers());
 202             }
 203           if (self->GetAbortExecute())
 204             {
 205             break;
 206             }
 207           }
 208
 209         // Set the integration step to be characteristic cell length
 210         // time IntegrationStepLength
 211         input->GetCell(pt1.cellId, cell);
 212         step = dir*self->GetIntegrationStepLength()
 213           * sqrt((double)cell->GetLength2())/pt1.speed;
 214
 215         // Calculate the next step using the integrator provided
 216         if (integrator->ComputeNextStep(pt1.x, pt1.v, xNext, 0, step, 0, err)
 217             != 0)
 218           {
 219           break;
 220           }
 221
 222         for(i=0; i<3; i++)
 223           {
 224           coords[i] = xNext[i];
 225           }
 226
 227         // Interpolate the velocity field at coords
 228         if ( !func->FunctionValues(coords, vel) )
 229           {
 230           break;
 231           }
 232
 233         for(i=0; i<3; i++)
 234           {
 235           pt2.v[i] = vel[i];
 236           }
 237
 238         for (i=0; i<3; i++)
 239           {
 240           pt2.x[i] = xNext[i];
 241           }
 242
 243         pt2.cellId = func->GetLastCellId();
 244         func->GetLastWeights(w);
 245         func->GetLastLocalCoordinates(pcoords);
 246         input->GetCell(pt2.cellId, cell);
 247
 248         if ( inScalars )
 249           {
 250           // Interpolate scalars
 251           inScalars->GetTuples(cell->PointIds, cellScalars);
 252           for (pt2.s=0.0, i=0; i < cell->GetNumberOfPoints(); i++)
 253             {
 254             pt2.s += cellScalars->GetComponent(i,0) * w[i];
 255             }
 256           }
 257
 258         pt2.speed = vtkMath::Norm(pt2.v);
 259
 260         d = sqrt((double)vtkMath::Distance2BetweenPoints(pt1.x,pt2.x));
 261         pt2.d = pt1.d + d;
 262         // If at stagnation region, stop the integration
 263         if ( d < VTK_EPSILON || (pt1.speed + pt2.speed) < VTK_EPSILON )
 264           {
 265           pt2.t = pt1.t;
 266           break;
 267           }
 268         pt2.t = pt1.t + (2.0 * d / (pt1.speed + pt2.speed));
 269
 270         if (self->GetVorticity() && inVectors)
 271           {
 272           // compute vorticity
 273           inVectors->GetTuples(cell->PointIds, cellVectors);
 274           cellVel = cellVectors->GetPointer(0);
 275           cell->Derivatives(0, pcoords, cellVel, 3, derivs);
 276           vort[0] = derivs[7] - derivs[5];
 277           vort[1] = derivs[2] - derivs[6];
 278           vort[2] = derivs[3] - derivs[1];
 279           // rotation
 280           pt2.omega = vtkMath::Dot(vort, pt2.v);
 281           pt2.omega /= pt2.speed;
 282           pt2.theta += (pt1.omega+pt2.omega)/2 * (pt2.t - pt1.t);
 283           }
 284
 285
 286         // Store only points which have a point to be displayed
 287         // between them
 288         if (tOffset >= pt1.t && tOffset <= pt2.t)
 289           {
 290           // Do not store if same as the last point.
 291           // To avoid storing some points twice.
 292           if ( !sNext || sNext->x[0] != pt1.x[0] || sNext->x[1] != pt1.x[1]
 293                || sNext->x[2] != pt1.x[2] )
 294             {
 295             idxNext = streamer->InsertNextStreamPoint();
 296             sNext = streamer->GetStreamPoint(idxNext);
 297             *sNext = pt1;
 298             }
 299           idxNext = streamer->InsertNextStreamPoint();
 300           sNext = streamer->GetStreamPoint(idxNext);
 301           *sNext = pt2;
 302           }
 303         if (tOffset < pt2.t)
 304           {
 305           tOffset += ((int)(( pt2.t - tOffset) / savePointInterval) + 1)
 306             * savePointInterval;
 307           }
 308         pt1 = pt2;
 309
 310         }
 311       // Store the last point anyway.
 312       if ( !sNext || sNext->x[0] != pt2.x[0] || sNext->x[1] != pt2.x[1]
 313            || sNext->x[2] != pt2.x[2] )
 314         {
 315         idxNext = streamer->InsertNextStreamPoint();
 316         sNext = streamer->GetStreamPoint(idxNext);
 317         *sNext = pt2;
 318         }
 319       // Clear the last cell to avoid starting a search from
 320       // the last point in the streamline
 321       func->ClearLastCellId();
 322       }
 323     }
 324
 325   integrator->Delete();
 326   func->Delete();
 327
 328   cell->Delete();
 329   cellVectors->Delete();
 330   if (cellScalars)
 331     {
 332     cellScalars->Delete();
 333     }
 334   delete[] w;
 335
 336   return VTK_THREAD_RETURN_VALUE;
 337 }
 338
 339 void VISU_Streamer::Integrate(vtkDataSet *input, vtkDataSet *source)
 340 {
 341   vtkPointData *pd   = input->GetPointData();
 342   vtkDataArray *inScalars;
 343   vtkDataArray *inVectors;
 344   vtkIdType numSourcePts, idx, idxNext;
 345   vtkStreamer::StreamPoint *sNext, *sPtr;
 346   vtkIdType ptId, i;
 347   int j, offset;
 348   vtkCell *cell;
 349   double v[3], *cellVel, derivs[9], xNext[3], vort[3];
 350   double tol2;
 351   double *w = new double[input->GetMaxCellSize()];
 352   vtkDoubleArray *cellVectors;
 353   vtkDataArray *cellScalars=0;
 354
 355   vtkDebugMacro(<<"Generating streamers");
 356   this->NumberOfStreamers = 0;
 357
 358   // reexecuting - delete old stuff
 359   if( this->Streamers )
 360     {
 361     delete [] this->Streamers;
 362     }
 363   this->Streamers = NULL;
 364
 365   if ( ! (inVectors=pd->GetVectors()) )
 366     {
 367     delete [] w;
 368     vtkErrorMacro(<<"No vector data defined!");
 369     return;
 370     }
 371
 372   cellVectors = vtkDoubleArray::New();
 373   cellVectors->SetNumberOfComponents(3);
 374   cellVectors->Allocate(3*VTK_CELL_SIZE);
 375
 376   inScalars = pd->GetScalars();
 377
 378   if (inScalars)
 379     {
 380     cellScalars = inScalars->NewInstance();
 381     cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
 382     cellScalars->Allocate(cellScalars->GetNumberOfComponents()*VTK_CELL_SIZE);
 383     }
 384
 385   tol2 = input->GetLength()/1000;
 386   tol2 = tol2*tol2;
 387
 388   //
 389   // Create starting points
 390   //
 391   this->NumberOfStreamers = numSourcePts = offset = 1;
 392   if ( source )
 393     {
 394     this->NumberOfStreamers = numSourcePts = source->GetNumberOfPoints();
 395     }
 396
 397   if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
 398     {
 399     offset = 2;
 400     this->NumberOfStreamers *= 2;
 401     }
 402
 403   this->Streamers = new vtkStreamer::StreamArray[this->NumberOfStreamers];
 404
 405   if ( this->StartFrom == VTK_START_FROM_POSITION && !source )
 406     {
 407     idx = this->Streamers[0].InsertNextStreamPoint();
 408     sPtr = this->Streamers[0].GetStreamPoint(idx);
 409     sPtr->subId = 0;
 410     for (i=0; i<3; i++)
 411       {
 412       sPtr->x[i] = this->StartPosition[i];
 413       }
 414     sPtr->cellId = input->FindCell(this->StartPosition, NULL, -1, 0.0,
 415                                    sPtr->subId, sPtr->p, w);
 416     }
 417
 418   else if ( this->StartFrom == VTK_START_FROM_LOCATION && !source )
 419     {
 420     idx = this->Streamers[0].InsertNextStreamPoint();
 421     sPtr = this->Streamers[0].GetStreamPoint(idx);
 422     sPtr->subId = 0;
 423     cell =  input->GetCell(sPtr->cellId);
 424     cell->EvaluateLocation(sPtr->subId, sPtr->p, sPtr->x, w);
 425     }
 426
 427   else //VTK_START_FROM_SOURCE
 428     {
 429     for (ptId=0; ptId < numSourcePts; ptId++)
 430       {
 431       idx = this->Streamers[offset*ptId].InsertNextStreamPoint();
 432       sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
 433       sPtr->subId = 0;
 434       source->GetPoint(ptId,sPtr->x);
 435       sPtr->cellId = input->FindCell(sPtr->x, NULL, -1, tol2,
 436                                      sPtr->subId, sPtr->p, w);
 437       }
 438     }
 439
 440   // Finish initializing each streamer
 441   //
 442   for (idx=0, ptId=0; ptId < numSourcePts; ptId++)
 443     {
 444     this->Streamers[offset*ptId].Direction = 1.0;
 445     sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
 446     sPtr->d = 0.0;
 447     sPtr->t = 0.0;
 448     sPtr->s = 0.0;
 449     sPtr->theta = 0.0;
 450     sPtr->omega = 0.0;
 451
 452     if ( sPtr->cellId >= 0 ) //starting point in dataset
 453       {
 454       cell = input->GetCell(sPtr->cellId);
 455       cell->EvaluateLocation(sPtr->subId, sPtr->p, xNext, w);
 456
 457       inVectors->GetTuples(cell->PointIds, cellVectors);
 458       sPtr->v[0]  = sPtr->v[1] = sPtr->v[2] = 0.0;
 459       for (i=0; i < cell->GetNumberOfPoints(); i++)
 460         {
 461         cellVectors->GetTuple(i, v);
 462         for (j=0; j<3; j++)
 463           {
 464           sPtr->v[j] += v[j] * w[i];
 465           }
 466         }
 467
 468       sPtr->speed = vtkMath::Norm(sPtr->v);
 469
 470       if (this->GetVorticity() && inVectors)
 471         {
 472           // compute vorticity
 473         inVectors->GetTuples(cell->PointIds, cellVectors);
 474         cellVel = cellVectors->GetPointer(0);
 475         cell->Derivatives(0, sPtr->p, cellVel, 3, derivs);
 476         vort[0] = derivs[7] - derivs[5];
 477         vort[1] = derivs[2] - derivs[6];
 478         vort[2] = derivs[3] - derivs[1];
 479         // rotation
 480         sPtr->omega = vtkMath::Dot(vort, sPtr->v);
 481         sPtr->omega /= sPtr->speed;
 482         sPtr->theta = 0;
 483         }
 484
 485       if ( inScalars )
 486         {
 487         inScalars->GetTuples(cell->PointIds, cellScalars);
 488         for (sPtr->s=0, i=0; i < cell->GetNumberOfPoints(); i++)
 489           {
 490           sPtr->s += cellScalars->GetComponent(i,0) * w[i];
 491           }
 492         }
 493       }
 494     else
 495       {
 496       for (j=0; j<3; j++)
 497         {
 498         sPtr->p[j] = 0.0;
 499         sPtr->v[j] = 0.0;
 500         }
 501       sPtr->speed = 0;
 502       }
 503
 504     if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
 505       {
 506       this->Streamers[offset*ptId+1].Direction = -1.0;
 507       idxNext = this->Streamers[offset*ptId+1].InsertNextStreamPoint();
 508       sNext = this->Streamers[offset*ptId+1].GetStreamPoint(idxNext);
 509       sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
 510       *sNext = *sPtr;
 511       }
 512     else if ( this->IntegrationDirection == VTK_INTEGRATE_BACKWARD )
 513       {
 514       this->Streamers[offset*ptId].Direction = -1.0;
 515       }
 516
 517
 518     } //for each streamer
 519
 520   // Some data access methods must be called once from a single thread before they
 521   // can safely be used. Call those now
 522   vtkGenericCell *gcell = vtkGenericCell::New();
 523   input->GetCell(0,gcell);
 524   gcell->Delete();
 525
 526   // Set up and execute the thread
 527   this->Threader->SetNumberOfThreads( this->NumberOfThreads );
 528   VISU_StreamerThreadStruct str;
 529   str.Filter = this;
 530   str.Input = input;
 531   str.Source = source;
 532   this->Threader->SetSingleMethod( VISU_Streamer::ThreadedIntegrate, &str );
 533   this->Threader->SingleMethodExecute();
 534
 535   //
 536   // Now create appropriate representation
 537   //
 538   if ( this->OrientationScalars && !this->SpeedScalars)
 539     {
 540     for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
 541       {
 542       for ( sPtr=this->Streamers[ptId].GetStreamPoint(0), i=0;
 543             i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
 544             i++, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
 545         {
 546         sPtr->s = sPtr->theta;
 547         }
 548       }
 549     }
 550
 551   if ( this->SpeedScalars )
 552     {
 553     for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
 554       {
 555       for ( sPtr=this->Streamers[ptId].GetStreamPoint(0), i=0;
 556             i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
 557             i++, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
 558         {
 559         sPtr->s = sPtr->speed;
 560         }
 561       }
 562     }
 563   delete [] w;
 564   cellVectors->Delete();
 565   if (cellScalars)
 566     {
 567     cellScalars->Delete();
 568     }
 569 }