1 // Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "VISU_Streamer.hxx"
23 #include "vtkDataSet.h"
24 #include "vtkDoubleArray.h"
25 #include "vtkExecutive.h"
26 #include "vtkGenericCell.h"
27 #include "vtkInformation.h"
28 #include "vtkInterpolatedVelocityField.h"
30 #include "vtkMultiThreader.h"
31 #include "vtkObjectFactory.h"
32 #include "vtkPointData.h"
33 #include "vtkRungeKutta2.h"
35 vtkCxxRevisionMacro(VISU_Streamer, "$Revision$");
37 #define VTK_START_FROM_POSITION 0
38 #define VTK_START_FROM_LOCATION 1
40 struct VISU_StreamerThreadStruct
42 VISU_Streamer *Filter;
47 vtkStreamer::StreamArray::StreamArray()
50 this->Array = new vtkStreamer::StreamPoint[1000];
53 this->Direction = VTK_INTEGRATE_FORWARD;
56 vtkStreamer::StreamPoint *vtkStreamer::StreamArray::Resize(vtkIdType sz)
58 vtkStreamer::StreamPoint *newArray;
63 newSize = this->Size +
64 this->Extend*(((sz-this->Size)/this->Extend)+1);
71 newArray = new vtkStreamer::StreamPoint[newSize];
73 memcpy(newArray, this->Array,
74 (sz < this->Size ? sz : this->Size) * sizeof(vtkStreamer::StreamPoint));
77 delete [] this->Array;
78 this->Array = newArray;
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()
90 VISU_Streamer::~VISU_Streamer()
94 static const double VTK_EPSILON=1E-7;
96 VTK_THREAD_RETURN_TYPE VISU_Streamer::ThreadedIntegrate( void *arg )
99 VISU_StreamerThreadStruct *str;
102 vtkStreamer::StreamArray *streamer;
103 vtkStreamer::StreamPoint *sNext = 0, *sPtr;
104 vtkStreamer::StreamPoint pt1, pt2;
106 vtkIdType idxNext, ptId;
108 double xNext[3], vel[3];
111 double *w, pcoords[3];
114 vtkGenericCell *cell;
116 vtkDataArray *inScalars;
117 vtkDataArray *inVectors;
118 vtkDoubleArray *cellVectors;
119 vtkDataArray *cellScalars=0;
120 double tOffset, vort[3];
124 thread_id = ((vtkMultiThreader::ThreadInfo *)(arg))->ThreadID;
125 thread_count = ((vtkMultiThreader::ThreadInfo *)(arg))->NumberOfThreads;
126 str = (VISU_StreamerThreadStruct*)(((vtkMultiThreader::ThreadInfo *)(arg))->UserData);
130 pd = input->GetPointData();
131 inScalars = pd->GetScalars();
132 inVectors = pd->GetVectors();
134 cell = vtkGenericCell::New();
135 cellVectors = vtkDoubleArray::New();
136 cellVectors->SetNumberOfComponents(3);
137 cellVectors->Allocate(3*VTK_CELL_SIZE);
140 cellScalars = inScalars->NewInstance();
141 cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
142 cellScalars->Allocate(inScalars->GetNumberOfComponents()*VTK_CELL_SIZE);
145 w = new double[input->GetMaxCellSize()];
147 // Set the function set to be integrated
148 vtkInterpolatedVelocityField* func = vtkInterpolatedVelocityField::New();
149 func->AddDataSet(input);
151 if (self->GetIntegrator() == 0)
153 vtkGenericWarningMacro("No integrator is specified.");
154 return VTK_THREAD_RETURN_VALUE;
157 // Create a new integrator, the type is the same as Integrator
158 vtkInitialValueProblemSolver* integrator =
159 self->GetIntegrator()->NewInstance();
160 integrator->SetFunctionSet(func);
162 // Used to avoid calling these function many times during
164 double termspeed = self->GetTerminalSpeed();
165 double maxtime = self->GetMaximumPropagationTime();
166 double savePointInterval = self->GetSavePointInterval();
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++)
172 if ( ptId % thread_count == thread_id )
175 streamer = self->GetStreamers() + ptId;
176 sPtr = streamer->GetStreamPoint(0);
177 if ( sPtr->cellId < 0 )
181 // Set the last cell id in the vtkInterpolatedVelocityField
182 // object to speed up FindCell calls
183 func->SetLastCellId(sPtr->cellId);
185 dir = streamer->Direction;
187 // Copy the first point
192 //integrate until time has been exceeded
193 while ( pt1.cellId >= 0 && pt1.speed > termspeed && pt1.t < maxtime )
196 if ( counter++ % 1000 == 0 )
200 self->UpdateProgress((double)ptId/self->GetNumberOfStreamers()
201 +pt1.t/maxtime/self->GetNumberOfStreamers());
203 if (self->GetAbortExecute())
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;
215 // Calculate the next step using the integrator provided
216 if (integrator->ComputeNextStep(pt1.x, pt1.v, xNext, 0, step, 0, err)
224 coords[i] = xNext[i];
227 // Interpolate the velocity field at coords
228 if ( !func->FunctionValues(coords, vel) )
243 pt2.cellId = func->GetLastCellId();
244 func->GetLastWeights(w);
245 func->GetLastLocalCoordinates(pcoords);
246 input->GetCell(pt2.cellId, cell);
250 // Interpolate scalars
251 inScalars->GetTuples(cell->PointIds, cellScalars);
252 for (pt2.s=0.0, i=0; i < cell->GetNumberOfPoints(); i++)
254 pt2.s += cellScalars->GetComponent(i,0) * w[i];
258 pt2.speed = vtkMath::Norm(pt2.v);
260 d = sqrt((double)vtkMath::Distance2BetweenPoints(pt1.x,pt2.x));
262 // If at stagnation region, stop the integration
263 if ( d < VTK_EPSILON || (pt1.speed + pt2.speed) < VTK_EPSILON )
268 pt2.t = pt1.t + (2.0 * d / (pt1.speed + pt2.speed));
270 if (self->GetVorticity() && inVectors)
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];
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);
286 // Store only points which have a point to be displayed
288 if (tOffset >= pt1.t && tOffset <= pt2.t)
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] )
295 idxNext = streamer->InsertNextStreamPoint();
296 sNext = streamer->GetStreamPoint(idxNext);
299 idxNext = streamer->InsertNextStreamPoint();
300 sNext = streamer->GetStreamPoint(idxNext);
305 tOffset += ((int)(( pt2.t - tOffset) / savePointInterval) + 1)
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] )
315 idxNext = streamer->InsertNextStreamPoint();
316 sNext = streamer->GetStreamPoint(idxNext);
319 // Clear the last cell to avoid starting a search from
320 // the last point in the streamline
321 func->ClearLastCellId();
325 integrator->Delete();
329 cellVectors->Delete();
332 cellScalars->Delete();
336 return VTK_THREAD_RETURN_VALUE;
339 void VISU_Streamer::Integrate(vtkDataSet *input, vtkDataSet *source)
341 vtkPointData *pd = input->GetPointData();
342 vtkDataArray *inScalars;
343 vtkDataArray *inVectors;
344 vtkIdType numSourcePts, idx, idxNext;
345 vtkStreamer::StreamPoint *sNext, *sPtr;
349 double v[3], *cellVel, derivs[9], xNext[3], vort[3];
351 double *w = new double[input->GetMaxCellSize()];
352 vtkDoubleArray *cellVectors;
353 vtkDataArray *cellScalars=0;
355 vtkDebugMacro(<<"Generating streamers");
356 this->NumberOfStreamers = 0;
358 // reexecuting - delete old stuff
359 if( this->Streamers )
361 delete [] this->Streamers;
363 this->Streamers = NULL;
365 if ( ! (inVectors=pd->GetVectors()) )
368 vtkErrorMacro(<<"No vector data defined!");
372 cellVectors = vtkDoubleArray::New();
373 cellVectors->SetNumberOfComponents(3);
374 cellVectors->Allocate(3*VTK_CELL_SIZE);
376 inScalars = pd->GetScalars();
380 cellScalars = inScalars->NewInstance();
381 cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
382 cellScalars->Allocate(cellScalars->GetNumberOfComponents()*VTK_CELL_SIZE);
385 tol2 = input->GetLength()/1000;
389 // Create starting points
391 this->NumberOfStreamers = numSourcePts = offset = 1;
394 this->NumberOfStreamers = numSourcePts = source->GetNumberOfPoints();
397 if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
400 this->NumberOfStreamers *= 2;
403 this->Streamers = new vtkStreamer::StreamArray[this->NumberOfStreamers];
405 if ( this->StartFrom == VTK_START_FROM_POSITION && !source )
407 idx = this->Streamers[0].InsertNextStreamPoint();
408 sPtr = this->Streamers[0].GetStreamPoint(idx);
412 sPtr->x[i] = this->StartPosition[i];
414 sPtr->cellId = input->FindCell(this->StartPosition, NULL, -1, 0.0,
415 sPtr->subId, sPtr->p, w);
418 else if ( this->StartFrom == VTK_START_FROM_LOCATION && !source )
420 idx = this->Streamers[0].InsertNextStreamPoint();
421 sPtr = this->Streamers[0].GetStreamPoint(idx);
423 cell = input->GetCell(sPtr->cellId);
424 cell->EvaluateLocation(sPtr->subId, sPtr->p, sPtr->x, w);
427 else //VTK_START_FROM_SOURCE
429 for (ptId=0; ptId < numSourcePts; ptId++)
431 idx = this->Streamers[offset*ptId].InsertNextStreamPoint();
432 sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
434 source->GetPoint(ptId,sPtr->x);
435 sPtr->cellId = input->FindCell(sPtr->x, NULL, -1, tol2,
436 sPtr->subId, sPtr->p, w);
440 // Finish initializing each streamer
442 for (idx=0, ptId=0; ptId < numSourcePts; ptId++)
444 this->Streamers[offset*ptId].Direction = 1.0;
445 sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
452 if ( sPtr->cellId >= 0 ) //starting point in dataset
454 cell = input->GetCell(sPtr->cellId);
455 cell->EvaluateLocation(sPtr->subId, sPtr->p, xNext, w);
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++)
461 cellVectors->GetTuple(i, v);
464 sPtr->v[j] += v[j] * w[i];
468 sPtr->speed = vtkMath::Norm(sPtr->v);
470 if (this->GetVorticity() && inVectors)
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];
480 sPtr->omega = vtkMath::Dot(vort, sPtr->v);
481 sPtr->omega /= sPtr->speed;
487 inScalars->GetTuples(cell->PointIds, cellScalars);
488 for (sPtr->s=0, i=0; i < cell->GetNumberOfPoints(); i++)
490 sPtr->s += cellScalars->GetComponent(i,0) * w[i];
504 if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
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);
512 else if ( this->IntegrationDirection == VTK_INTEGRATE_BACKWARD )
514 this->Streamers[offset*ptId].Direction = -1.0;
518 } //for each streamer
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);
526 // Set up and execute the thread
527 this->Threader->SetNumberOfThreads( this->NumberOfThreads );
528 VISU_StreamerThreadStruct str;
532 this->Threader->SetSingleMethod( VISU_Streamer::ThreadedIntegrate, &str );
533 this->Threader->SingleMethodExecute();
536 // Now create appropriate representation
538 if ( this->OrientationScalars && !this->SpeedScalars)
540 for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
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) )
546 sPtr->s = sPtr->theta;
551 if ( this->SpeedScalars )
553 for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
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) )
559 sPtr->s = sPtr->speed;
564 cellVectors->Delete();
567 cellScalars->Delete();