VISU_IsoSurfacesPL.hxx \
VISU_DeformedShapePL.hxx \
VISU_VectorsPL.hxx \
+ VISU_Streamer.hxx \
+ VISU_StreamLine.hxx \
VISU_StreamLinesPL.hxx \
VISU_LookupTable.hxx \
VISU_ScalarBarActor.hxx \
VISU_IsoSurfacesPL.cxx \
VISU_DeformedShapePL.cxx \
VISU_VectorsPL.cxx \
+ VISU_Streamer.cxx \
+ VISU_StreamLine.cxx \
VISU_StreamLinesPL.cxx \
VISU_LookupTable.cxx \
VISU_ScalarBarActor.cxx \
//----------------------------------------------------------------------------
size_t
VISU_PipeLine
-::GetAvailableMemory(size_t theSize,
- size_t theMinSize)
+::GetAvailableMemory(double theSize,
+ double theMinSize)
{
// Finds acceptable memory size by half-deflection methods
static double EPSILON = 2 * 1024;
static
size_t
- GetAvailableMemory(size_t theSize,
- size_t theMinSize = 1024*1024);
+ GetAvailableMemory(double theSize,
+ double theMinSize = 1024*1024);
protected:
//----------------------------------------------------------------------------
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#include "VISU_StreamLine.hxx"
+
+#include "vtkCellArray.h"
+#include "vtkDataSet.h"
+#include "vtkFloatArray.h"
+#include "vtkInformation.h"
+#include "vtkInformationVector.h"
+#include "vtkMath.h"
+#include "vtkObjectFactory.h"
+#include "vtkPointData.h"
+#include "vtkPolyData.h"
+#include "vtkPolyLine.h"
+
+vtkCxxRevisionMacro(VISU_StreamLine, "$Revision$");
+vtkStandardNewMacro(VISU_StreamLine);
+
+// Construct object with step size set to 1.0.
+VISU_StreamLine::VISU_StreamLine()
+{
+ this->StepLength = 1.0;
+ this->NumberOfStreamers = 0;
+}
+
+int VISU_StreamLine::RequestData(
+ vtkInformation *,
+ vtkInformationVector **inputVector,
+ vtkInformationVector *outputVector)
+{
+ vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
+ vtkInformation *outInfo = outputVector->GetInformationObject(0);
+ vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
+
+ vtkDataSet *input = vtkDataSet::SafeDownCast(
+ inInfo->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData *output = vtkPolyData::SafeDownCast(
+ outInfo->Get(vtkDataObject::DATA_OBJECT()));
+ vtkDataSet *source = 0;
+ if (sourceInfo)
+ {
+ source = vtkDataSet::SafeDownCast(
+ sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
+ }
+
+ vtkStreamer::StreamPoint *sPrev, *sPtr;
+ vtkPoints *newPts;
+ vtkFloatArray *newVectors;
+ vtkFloatArray *newScalars=NULL;
+ vtkCellArray *newLines;
+ vtkIdType ptId, i, id;
+ int j;
+ vtkIdList *pts;
+ double tOffset, x[3], v[3], s, r;
+ double theta;
+ vtkPolyLine* lineNormalGenerator = NULL;
+ vtkFloatArray* normals = NULL;
+ vtkFloatArray* rotation = 0;
+
+ this->SavePointInterval = this->StepLength;
+ this->VISU_Streamer::Integrate(input, source);
+ if ( this->NumberOfStreamers <= 0 ) {return 1;}
+
+ pts = vtkIdList::New();
+ pts->Allocate(2500);
+
+ //
+ // Convert streamer into lines. Lines may be dashed.
+ //
+ newPts = vtkPoints::New();
+ newPts->Allocate(1000);
+ newVectors = vtkFloatArray::New();
+ newVectors->SetNumberOfComponents(3);
+ newVectors->Allocate(3000);
+ if ( this->Vorticity )
+ {
+ lineNormalGenerator = vtkPolyLine::New();
+ normals = vtkFloatArray::New();
+ normals->SetNumberOfComponents(3);
+ normals->Allocate(3000);
+ rotation = vtkFloatArray::New();
+ rotation->SetNumberOfComponents(1);
+ rotation->Allocate(1000);
+ rotation->SetName("Thetas");
+ output->GetPointData()->AddArray(rotation);
+ }
+
+ if ( input->GetPointData()->GetScalars() || this->SpeedScalars
+ || this->OrientationScalars)
+ {
+ newScalars = vtkFloatArray::New();
+ newScalars->Allocate(1000);
+ }
+ newLines = vtkCellArray::New();
+ newLines->Allocate(newLines->EstimateSize(2*this->NumberOfStreamers,
+ VTK_CELL_SIZE));
+ //
+ // Loop over all streamers generating points
+ //
+ for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
+ {
+ if ( this->Streamers[ptId].GetNumberOfPoints() < 2 )
+ {
+ continue;
+ }
+ sPrev = this->Streamers[ptId].GetStreamPoint(0);
+ sPtr = this->Streamers[ptId].GetStreamPoint(1);
+
+ if ( this->Streamers[ptId].GetNumberOfPoints() == 2 && sPtr->cellId >= 0 )
+ {
+ continue;
+ }
+
+ tOffset = sPrev->t;
+
+ for ( i=1;
+ i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
+ i++, sPrev=sPtr, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
+ {
+ //
+ // Create points for line
+ //
+ while ( tOffset >= sPrev->t && tOffset < sPtr->t )
+ {
+ r = (tOffset - sPrev->t) / (sPtr->t - sPrev->t);
+
+ for (j=0; j<3; j++)
+ {
+ x[j] = sPrev->x[j] + r * (sPtr->x[j] - sPrev->x[j]);
+ v[j] = sPrev->v[j] + r * (sPtr->v[j] - sPrev->v[j]);
+ }
+
+ // add point to line
+ id = newPts->InsertNextPoint(x);
+ pts->InsertNextId(id);
+ newVectors->InsertTuple(id,v);
+
+ if ( newScalars )
+ {
+ s = sPrev->s + r * (sPtr->s - sPrev->s);
+ newScalars->InsertTuple(id,&s);
+ }
+
+ if ( this->Vorticity )
+ {
+ // Store the rotation values. Used after all the streamlines
+ // are generated.
+ theta = sPrev->theta + r * (sPtr->theta - sPrev->theta);
+ rotation->InsertTuple(id, &theta);
+ }
+
+ tOffset += this->StepLength;
+
+ } // while
+ } //for this streamer
+
+ if ( pts->GetNumberOfIds() > 1 )
+ {
+ newLines->InsertNextCell(pts);
+ pts->Reset();
+ }
+ } //for all streamers
+
+ vtkDebugMacro(<<"Created " << newPts->GetNumberOfPoints() << " points, "
+ << newLines->GetNumberOfCells() << " lines");
+
+ if (this->Vorticity)
+ {
+ // Rotate the normal vectors with stream vorticity
+ vtkIdType nPts=0;
+ vtkIdType *linePts=0;
+ double normal[3], local1[3], local2[3], length, costheta, sintheta;
+
+ lineNormalGenerator->GenerateSlidingNormals(newPts,newLines,normals);
+
+ // Loop over all lines, from the above code we are know that each line
+ // will have at least two points and that no points will be shared
+ // between lines. It is important to loop over the points used by the
+ // lines because newPts may actually contain points that are not used by
+ // any lines. The normals are only calculated for points that are used
+ // in lines so referencing normals for all points can lead to UMRs
+ for (newLines->InitTraversal(); newLines->GetNextCell(nPts,linePts); )
+ {
+ for(i=0; i<nPts; i++)
+ {
+ normals->GetTuple(linePts[i], normal);
+ newVectors->GetTuple(linePts[i], v);
+ // obtain two unit orthogonal vectors on the plane perpendicular to
+ // the streamline
+ for(j=0; j<3; j++)
+ {
+ local1[j] = normal[j];
+ }
+ length = vtkMath::Normalize(local1);
+ vtkMath::Cross(local1, v, local2);
+ vtkMath::Normalize(local2);
+ // Rotate the normal with theta
+ rotation->GetTuple(linePts[i], &theta);
+ costheta = cos(theta);
+ sintheta = sin(theta);
+ for(j=0; j<3; j++)
+ {
+ normal[j] = length* (costheta*local1[j] + sintheta*local2[j]);
+ }
+ normals->SetTuple(linePts[i], normal);
+ }
+ }
+ output->GetPointData()->SetNormals(normals);
+ normals->Delete();
+ lineNormalGenerator->Delete();
+ rotation->Delete();
+ }
+
+ output->SetPoints(newPts);
+ newPts->Delete();
+
+ output->GetPointData()->SetVectors(newVectors);
+ newVectors->Delete();
+
+ if ( newScalars )
+ {
+ int idx = output->GetPointData()->AddArray(newScalars);
+ output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
+ newScalars->Delete();
+ }
+
+ pts->Delete();
+ output->SetLines(newLines);
+ newLines->Delete();
+
+ // Delete the streamers since they are no longer needed
+ delete[] this->Streamers;
+ this->Streamers = 0;
+ this->NumberOfStreamers = 0;
+
+ output->Squeeze();
+
+ return 1;
+}
+
+void VISU_StreamLine::PrintSelf(ostream& os, vtkIndent indent)
+{
+ this->Superclass::PrintSelf(os,indent);
+
+ os << indent << "Step Length: " << this->StepLength << "\n";
+
+}
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#ifndef VISU_StreamLine_HeaderFile
+#define VISU_StreamLine_HeaderFile
+
+#include "VISUPipeline.hxx"
+#include "VISU_Streamer.hxx"
+
+class VISU_PIPELINE_EXPORT VISU_StreamLine : public VISU_Streamer
+{
+public:
+ vtkTypeRevisionMacro(VISU_StreamLine,VISU_Streamer);
+ void PrintSelf(ostream& os, vtkIndent indent);
+
+ // Description:
+ // Construct object with step size set to 1.0.
+ static VISU_StreamLine *New();
+
+ // Description:
+ // Specify the length of a line segment. The length is expressed in terms of
+ // elapsed time. Smaller values result in smoother appearing streamlines, but
+ // greater numbers of line primitives.
+ vtkSetClampMacro(StepLength,double,0.000001,VTK_DOUBLE_MAX);
+ vtkGetMacro(StepLength,double);
+
+protected:
+ VISU_StreamLine();
+ ~VISU_StreamLine() {};
+
+ // Convert streamer array into vtkPolyData
+ virtual int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *);
+
+ // the length of line primitives
+ double StepLength;
+
+private:
+ VISU_StreamLine(const VISU_StreamLine&); // Not implemented.
+ void operator=(const VISU_StreamLine&); // Not implemented.
+};
+
+#endif
//#include "VISU_UsedPointsFilter.hxx"
#include "VISU_MaskPointsFilter.hxx"
#include "VISU_PipeLineUtils.hxx"
+#include "VISU_StreamLine.hxx"
#include "VTKViewer_GeometryFilter.h"
#include <vtkCell.h>
#include <vtkDataSet.h>
-#include <vtkStreamLine.h>
#ifdef _DEBUG_
static int MYDEBUG = 0;
SetIsShrinkable(false);
SetIsFeatureEdgesAllowed(false);
- myStream = vtkStreamLine::New();
+ myStream = VISU_StreamLine::New();
myCenters = vtkCellCenters::New();
myGeomFilter = VTKViewer_GeometryFilter::New();
myPointsFilter = VISU_MaskPointsFilter::New();
//----------------------------------------------------------------------------
-size_t
+vtkFloatingPointType
VISU_StreamLinesPL
::GetNecasseryMemorySize(vtkIdType theNbOfPoints,
vtkFloatingPointType theStepLength,
vtkFloatingPointType anAssignedDataSize = aCellsSize*4.0*sizeof(vtkFloatingPointType);
vtkFloatingPointType anOutputDataSetSize = aMeshSize + anAssignedDataSize;
- size_t aResult = size_t(aStreamArraySize*aNbCells + anOutputDataSetSize);
+ vtkFloatingPointType aResult = aStreamArraySize*aNbCells + anOutputDataSetSize;
return aResult;
}
{
static vtkFloatingPointType aPercentsDecrease = 3.0, aStepLengthIncrease = 9.0;
vtkIdType aNbOfPoints = theDataSet->GetNumberOfPoints();
- size_t aSize = GetNecasseryMemorySize(aNbOfPoints,theStepLength,thePropogationTime,thePercents);
+ vtkFloatingPointType aSize = GetNecasseryMemorySize(aNbOfPoints,theStepLength,thePropogationTime,thePercents);
size_t anIsPossible = CheckAvailableMemory(aSize);
if(!anIsPossible){
vtkFloatingPointType aMaxStepLength = std::max(GetMaxStepLength(theDataSet),thePropogationTime);
return 0.0; // absolutely empty object
vtkFloatingPointType anStepLength = GetMaxIntegrationStep(theDataSet)/aCoeffOfIntStep;
- vtkFloatingPointType aBasePropTime = GetBasePropagationTime(theDataSet)/GetVelocityCoeff(theDataSet);
+ // 0020724: last division has been commented, seems to be a logical mistake (to discuss with APO)
+ vtkFloatingPointType aBasePropTime = GetBasePropagationTime(theDataSet); // /GetVelocityCoeff(theDataSet)
thePercents = 1.0;
vtkIdType aNbOfPoints = theDataSet->GetNumberOfPoints();
- size_t aSize = GetNecasseryMemorySize(aNbOfPoints,anStepLength,aBasePropTime,thePercents);
+ vtkFloatingPointType aSize = GetNecasseryMemorySize(aNbOfPoints,anStepLength,aBasePropTime,thePercents);
size_t aRealSize = GetAvailableMemory(aSize);
vtkFloatingPointType anAverageVolume = aVolume / aRealSize;
vtkFloatingPointType aStep = pow(double(anAverageVolume), double(1.0/double(degree)));
#include "VISUPipeline.hxx"
#include "VISU_DeformedShapePL.hxx"
-#include <vtkStreamLine.h>
+
+#include <vtkStreamer.h>
class vtkDataSet;
class vtkPointSet;
class vtkCellCenters;
class VTKViewer_GeometryFilter;
class VISU_MaskPointsFilter;
+class VISU_StreamLine;
//----------------------------------------------------------------------------
bool theIsCopyInput);
static
- size_t
+ vtkFloatingPointType
GetNecasseryMemorySize(vtkIdType theNbOfPoints,
vtkFloatingPointType theStepLength,
vtkFloatingPointType thePropogationTime,
CorrectStepLength(vtkFloatingPointType theStep,
vtkDataSet* theDataSet);
- vtkStreamLine* myStream;
+ VISU_StreamLine* myStream;
vtkPointSet* mySource;
vtkCellCenters* myCenters;
VTKViewer_GeometryFilter *myGeomFilter;
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#include "VISU_Streamer.hxx"
+
+#include "vtkCell.h"
+#include "vtkDataSet.h"
+#include "vtkDoubleArray.h"
+#include "vtkExecutive.h"
+#include "vtkGenericCell.h"
+#include "vtkInformation.h"
+#include "vtkInterpolatedVelocityField.h"
+#include "vtkMath.h"
+#include "vtkMultiThreader.h"
+#include "vtkObjectFactory.h"
+#include "vtkPointData.h"
+#include "vtkRungeKutta2.h"
+
+vtkCxxRevisionMacro(VISU_Streamer, "$Revision$");
+
+#define VTK_START_FROM_POSITION 0
+#define VTK_START_FROM_LOCATION 1
+
+struct VISU_StreamerThreadStruct
+{
+ VISU_Streamer *Filter;
+ vtkDataSet *Input;
+ vtkDataSet *Source;
+};
+
+vtkStreamer::StreamArray::StreamArray()
+{
+ this->MaxId = -1;
+ this->Array = new vtkStreamer::StreamPoint[1000];
+ this->Size = 1000;
+ this->Extend = 5000;
+ this->Direction = VTK_INTEGRATE_FORWARD;
+}
+
+vtkStreamer::StreamPoint *vtkStreamer::StreamArray::Resize(vtkIdType sz)
+{
+ vtkStreamer::StreamPoint *newArray;
+ vtkIdType newSize;
+
+ if (sz >= this->Size)
+ {
+ newSize = this->Size +
+ this->Extend*(((sz-this->Size)/this->Extend)+1);
+ }
+ else
+ {
+ newSize = sz;
+ }
+
+ newArray = new vtkStreamer::StreamPoint[newSize];
+
+ memcpy(newArray, this->Array,
+ (sz < this->Size ? sz : this->Size) * sizeof(vtkStreamer::StreamPoint));
+
+ this->Size = newSize;
+ delete [] this->Array;
+ this->Array = newArray;
+
+ return this->Array;
+}
+
+// Construct object to start from position (0,0,0); integrate forward; terminal
+// speed 0.0; vorticity computation off; integrations step length 0.2; and
+// maximum propagation time 100.0.
+VISU_Streamer::VISU_Streamer()
+{
+}
+
+VISU_Streamer::~VISU_Streamer()
+{
+}
+
+static const double VTK_EPSILON=1E-7;
+
+VTK_THREAD_RETURN_TYPE VISU_Streamer::ThreadedIntegrate( void *arg )
+{
+ VISU_Streamer *self;
+ VISU_StreamerThreadStruct *str;
+ int thread_count;
+ int thread_id;
+ vtkStreamer::StreamArray *streamer;
+ vtkStreamer::StreamPoint *sNext = 0, *sPtr;
+ vtkStreamer::StreamPoint pt1, pt2;
+ int i;
+ vtkIdType idxNext, ptId;
+ double d, step, dir;
+ double xNext[3], vel[3];
+ double *cellVel;
+ double derivs[9];
+ double *w, pcoords[3];
+ double coords[4];
+ vtkDataSet *input;
+ vtkGenericCell *cell;
+ vtkPointData *pd;
+ vtkDataArray *inScalars;
+ vtkDataArray *inVectors;
+ vtkDoubleArray *cellVectors;
+ vtkDataArray *cellScalars=0;
+ double tOffset, vort[3];
+ double err;
+ int counter=0;
+
+ thread_id = ((vtkMultiThreader::ThreadInfo *)(arg))->ThreadID;
+ thread_count = ((vtkMultiThreader::ThreadInfo *)(arg))->NumberOfThreads;
+ str = (VISU_StreamerThreadStruct*)(((vtkMultiThreader::ThreadInfo *)(arg))->UserData);
+ self = str->Filter;
+
+ input = str->Input;
+ pd = input->GetPointData();
+ inScalars = pd->GetScalars();
+ inVectors = pd->GetVectors();
+
+ cell = vtkGenericCell::New();
+ cellVectors = vtkDoubleArray::New();
+ cellVectors->SetNumberOfComponents(3);
+ cellVectors->Allocate(3*VTK_CELL_SIZE);
+ if (inScalars)
+ {
+ cellScalars = inScalars->NewInstance();
+ cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
+ cellScalars->Allocate(inScalars->GetNumberOfComponents()*VTK_CELL_SIZE);
+ }
+
+ w = new double[input->GetMaxCellSize()];
+
+ // Set the function set to be integrated
+ vtkInterpolatedVelocityField* func = vtkInterpolatedVelocityField::New();
+ func->AddDataSet(input);
+
+ if (self->GetIntegrator() == 0)
+ {
+ vtkGenericWarningMacro("No integrator is specified.");
+ return VTK_THREAD_RETURN_VALUE;
+ }
+
+ // Create a new integrator, the type is the same as Integrator
+ vtkInitialValueProblemSolver* integrator =
+ self->GetIntegrator()->NewInstance();
+ integrator->SetFunctionSet(func);
+
+ // Used to avoid calling these function many times during
+ // the integration
+ double termspeed = self->GetTerminalSpeed();
+ double maxtime = self->GetMaximumPropagationTime();
+ double savePointInterval = self->GetSavePointInterval();
+
+ // For each streamer, integrate in appropriate direction
+ // Do only the streamers that this thread should handle.
+ for (ptId=0; ptId < self->GetNumberOfStreamers(); ptId++)
+ {
+ if ( ptId % thread_count == thread_id )
+ {
+ // Get starting step
+ streamer = self->GetStreamers() + ptId;
+ sPtr = streamer->GetStreamPoint(0);
+ if ( sPtr->cellId < 0 )
+ {
+ continue;
+ }
+ // Set the last cell id in the vtkInterpolatedVelocityField
+ // object to speed up FindCell calls
+ func->SetLastCellId(sPtr->cellId);
+
+ dir = streamer->Direction;
+
+ // Copy the first point
+ pt1 = *sPtr;
+ pt2 = *sPtr;
+ tOffset = pt1.t;
+
+ //integrate until time has been exceeded
+ while ( pt1.cellId >= 0 && pt1.speed > termspeed && pt1.t < maxtime )
+ {
+
+ if ( counter++ % 1000 == 0 )
+ {
+ if (!thread_id)
+ {
+ self->UpdateProgress((double)ptId/self->GetNumberOfStreamers()
+ +pt1.t/maxtime/self->GetNumberOfStreamers());
+ }
+ if (self->GetAbortExecute())
+ {
+ break;
+ }
+ }
+
+ // Set the integration step to be characteristic cell length
+ // time IntegrationStepLength
+ input->GetCell(pt1.cellId, cell);
+ step = dir*self->GetIntegrationStepLength()
+ * sqrt((double)cell->GetLength2())/pt1.speed;
+
+ // Calculate the next step using the integrator provided
+ if (integrator->ComputeNextStep(pt1.x, pt1.v, xNext, 0, step, 0, err)
+ != 0)
+ {
+ break;
+ }
+
+ for(i=0; i<3; i++)
+ {
+ coords[i] = xNext[i];
+ }
+
+ // Interpolate the velocity field at coords
+ if ( !func->FunctionValues(coords, vel) )
+ {
+ break;
+ }
+
+ for(i=0; i<3; i++)
+ {
+ pt2.v[i] = vel[i];
+ }
+
+ for (i=0; i<3; i++)
+ {
+ pt2.x[i] = xNext[i];
+ }
+
+ pt2.cellId = func->GetLastCellId();
+ func->GetLastWeights(w);
+ func->GetLastLocalCoordinates(pcoords);
+ input->GetCell(pt2.cellId, cell);
+
+ if ( inScalars )
+ {
+ // Interpolate scalars
+ inScalars->GetTuples(cell->PointIds, cellScalars);
+ for (pt2.s=0.0, i=0; i < cell->GetNumberOfPoints(); i++)
+ {
+ pt2.s += cellScalars->GetComponent(i,0) * w[i];
+ }
+ }
+
+ pt2.speed = vtkMath::Norm(pt2.v);
+
+ d = sqrt((double)vtkMath::Distance2BetweenPoints(pt1.x,pt2.x));
+ pt2.d = pt1.d + d;
+ // If at stagnation region, stop the integration
+ if ( d < VTK_EPSILON || (pt1.speed + pt2.speed) < VTK_EPSILON )
+ {
+ pt2.t = pt1.t;
+ break;
+ }
+ pt2.t = pt1.t + (2.0 * d / (pt1.speed + pt2.speed));
+
+ if (self->GetVorticity() && inVectors)
+ {
+ // compute vorticity
+ inVectors->GetTuples(cell->PointIds, cellVectors);
+ cellVel = cellVectors->GetPointer(0);
+ cell->Derivatives(0, pcoords, cellVel, 3, derivs);
+ vort[0] = derivs[7] - derivs[5];
+ vort[1] = derivs[2] - derivs[6];
+ vort[2] = derivs[3] - derivs[1];
+ // rotation
+ pt2.omega = vtkMath::Dot(vort, pt2.v);
+ pt2.omega /= pt2.speed;
+ pt2.theta += (pt1.omega+pt2.omega)/2 * (pt2.t - pt1.t);
+ }
+
+
+ // Store only points which have a point to be displayed
+ // between them
+ if (tOffset >= pt1.t && tOffset <= pt2.t)
+ {
+ // Do not store if same as the last point.
+ // To avoid storing some points twice.
+ if ( !sNext || sNext->x[0] != pt1.x[0] || sNext->x[1] != pt1.x[1]
+ || sNext->x[2] != pt1.x[2] )
+ {
+ idxNext = streamer->InsertNextStreamPoint();
+ sNext = streamer->GetStreamPoint(idxNext);
+ *sNext = pt1;
+ }
+ idxNext = streamer->InsertNextStreamPoint();
+ sNext = streamer->GetStreamPoint(idxNext);
+ *sNext = pt2;
+ }
+ if (tOffset < pt2.t)
+ {
+ tOffset += ((int)(( pt2.t - tOffset) / savePointInterval) + 1)
+ * savePointInterval;
+ }
+ pt1 = pt2;
+
+ }
+ // Store the last point anyway.
+ if ( !sNext || sNext->x[0] != pt2.x[0] || sNext->x[1] != pt2.x[1]
+ || sNext->x[2] != pt2.x[2] )
+ {
+ idxNext = streamer->InsertNextStreamPoint();
+ sNext = streamer->GetStreamPoint(idxNext);
+ *sNext = pt2;
+ }
+ // Clear the last cell to avoid starting a search from
+ // the last point in the streamline
+ func->ClearLastCellId();
+ }
+ }
+
+ integrator->Delete();
+ func->Delete();
+
+ cell->Delete();
+ cellVectors->Delete();
+ if (cellScalars)
+ {
+ cellScalars->Delete();
+ }
+ delete[] w;
+
+ return VTK_THREAD_RETURN_VALUE;
+}
+
+void VISU_Streamer::Integrate(vtkDataSet *input, vtkDataSet *source)
+{
+ vtkPointData *pd = input->GetPointData();
+ vtkDataArray *inScalars;
+ vtkDataArray *inVectors;
+ vtkIdType numSourcePts, idx, idxNext;
+ vtkStreamer::StreamPoint *sNext, *sPtr;
+ vtkIdType ptId, i;
+ int j, offset;
+ vtkCell *cell;
+ double v[3], *cellVel, derivs[9], xNext[3], vort[3];
+ double tol2;
+ double *w = new double[input->GetMaxCellSize()];
+ vtkDoubleArray *cellVectors;
+ vtkDataArray *cellScalars=0;
+
+ vtkDebugMacro(<<"Generating streamers");
+ this->NumberOfStreamers = 0;
+
+ // reexecuting - delete old stuff
+ if( this->Streamers )
+ {
+ delete [] this->Streamers;
+ }
+ this->Streamers = NULL;
+
+ if ( ! (inVectors=pd->GetVectors()) )
+ {
+ delete [] w;
+ vtkErrorMacro(<<"No vector data defined!");
+ return;
+ }
+
+ cellVectors = vtkDoubleArray::New();
+ cellVectors->SetNumberOfComponents(3);
+ cellVectors->Allocate(3*VTK_CELL_SIZE);
+
+ inScalars = pd->GetScalars();
+
+ if (inScalars)
+ {
+ cellScalars = inScalars->NewInstance();
+ cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
+ cellScalars->Allocate(cellScalars->GetNumberOfComponents()*VTK_CELL_SIZE);
+ }
+
+ tol2 = input->GetLength()/1000;
+ tol2 = tol2*tol2;
+
+ //
+ // Create starting points
+ //
+ this->NumberOfStreamers = numSourcePts = offset = 1;
+ if ( source )
+ {
+ this->NumberOfStreamers = numSourcePts = source->GetNumberOfPoints();
+ }
+
+ if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
+ {
+ offset = 2;
+ this->NumberOfStreamers *= 2;
+ }
+
+ this->Streamers = new vtkStreamer::StreamArray[this->NumberOfStreamers];
+
+ if ( this->StartFrom == VTK_START_FROM_POSITION && !source )
+ {
+ idx = this->Streamers[0].InsertNextStreamPoint();
+ sPtr = this->Streamers[0].GetStreamPoint(idx);
+ sPtr->subId = 0;
+ for (i=0; i<3; i++)
+ {
+ sPtr->x[i] = this->StartPosition[i];
+ }
+ sPtr->cellId = input->FindCell(this->StartPosition, NULL, -1, 0.0,
+ sPtr->subId, sPtr->p, w);
+ }
+
+ else if ( this->StartFrom == VTK_START_FROM_LOCATION && !source )
+ {
+ idx = this->Streamers[0].InsertNextStreamPoint();
+ sPtr = this->Streamers[0].GetStreamPoint(idx);
+ sPtr->subId = 0;
+ cell = input->GetCell(sPtr->cellId);
+ cell->EvaluateLocation(sPtr->subId, sPtr->p, sPtr->x, w);
+ }
+
+ else //VTK_START_FROM_SOURCE
+ {
+ for (ptId=0; ptId < numSourcePts; ptId++)
+ {
+ idx = this->Streamers[offset*ptId].InsertNextStreamPoint();
+ sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
+ sPtr->subId = 0;
+ source->GetPoint(ptId,sPtr->x);
+ sPtr->cellId = input->FindCell(sPtr->x, NULL, -1, tol2,
+ sPtr->subId, sPtr->p, w);
+ }
+ }
+
+ // Finish initializing each streamer
+ //
+ for (idx=0, ptId=0; ptId < numSourcePts; ptId++)
+ {
+ this->Streamers[offset*ptId].Direction = 1.0;
+ sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
+ sPtr->d = 0.0;
+ sPtr->t = 0.0;
+ sPtr->s = 0.0;
+ sPtr->theta = 0.0;
+ sPtr->omega = 0.0;
+
+ if ( sPtr->cellId >= 0 ) //starting point in dataset
+ {
+ cell = input->GetCell(sPtr->cellId);
+ cell->EvaluateLocation(sPtr->subId, sPtr->p, xNext, w);
+
+ inVectors->GetTuples(cell->PointIds, cellVectors);
+ sPtr->v[0] = sPtr->v[1] = sPtr->v[2] = 0.0;
+ for (i=0; i < cell->GetNumberOfPoints(); i++)
+ {
+ cellVectors->GetTuple(i, v);
+ for (j=0; j<3; j++)
+ {
+ sPtr->v[j] += v[j] * w[i];
+ }
+ }
+
+ sPtr->speed = vtkMath::Norm(sPtr->v);
+
+ if (this->GetVorticity() && inVectors)
+ {
+ // compute vorticity
+ inVectors->GetTuples(cell->PointIds, cellVectors);
+ cellVel = cellVectors->GetPointer(0);
+ cell->Derivatives(0, sPtr->p, cellVel, 3, derivs);
+ vort[0] = derivs[7] - derivs[5];
+ vort[1] = derivs[2] - derivs[6];
+ vort[2] = derivs[3] - derivs[1];
+ // rotation
+ sPtr->omega = vtkMath::Dot(vort, sPtr->v);
+ sPtr->omega /= sPtr->speed;
+ sPtr->theta = 0;
+ }
+
+ if ( inScalars )
+ {
+ inScalars->GetTuples(cell->PointIds, cellScalars);
+ for (sPtr->s=0, i=0; i < cell->GetNumberOfPoints(); i++)
+ {
+ sPtr->s += cellScalars->GetComponent(i,0) * w[i];
+ }
+ }
+ }
+ else
+ {
+ for (j=0; j<3; j++)
+ {
+ sPtr->p[j] = 0.0;
+ sPtr->v[j] = 0.0;
+ }
+ sPtr->speed = 0;
+ }
+
+ if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
+ {
+ this->Streamers[offset*ptId+1].Direction = -1.0;
+ idxNext = this->Streamers[offset*ptId+1].InsertNextStreamPoint();
+ sNext = this->Streamers[offset*ptId+1].GetStreamPoint(idxNext);
+ sPtr = this->Streamers[offset*ptId].GetStreamPoint(idx);
+ *sNext = *sPtr;
+ }
+ else if ( this->IntegrationDirection == VTK_INTEGRATE_BACKWARD )
+ {
+ this->Streamers[offset*ptId].Direction = -1.0;
+ }
+
+
+ } //for each streamer
+
+ // Some data access methods must be called once from a single thread before they
+ // can safely be used. Call those now
+ vtkGenericCell *gcell = vtkGenericCell::New();
+ input->GetCell(0,gcell);
+ gcell->Delete();
+
+ // Set up and execute the thread
+ this->Threader->SetNumberOfThreads( this->NumberOfThreads );
+ VISU_StreamerThreadStruct str;
+ str.Filter = this;
+ str.Input = input;
+ str.Source = source;
+ this->Threader->SetSingleMethod( VISU_Streamer::ThreadedIntegrate, &str );
+ this->Threader->SingleMethodExecute();
+
+ //
+ // Now create appropriate representation
+ //
+ if ( this->OrientationScalars && !this->SpeedScalars)
+ {
+ for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
+ {
+ for ( sPtr=this->Streamers[ptId].GetStreamPoint(0), i=0;
+ i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
+ i++, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
+ {
+ sPtr->s = sPtr->theta;
+ }
+ }
+ }
+
+ if ( this->SpeedScalars )
+ {
+ for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
+ {
+ for ( sPtr=this->Streamers[ptId].GetStreamPoint(0), i=0;
+ i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
+ i++, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
+ {
+ sPtr->s = sPtr->speed;
+ }
+ }
+ }
+ delete [] w;
+ cellVectors->Delete();
+ if (cellScalars)
+ {
+ cellScalars->Delete();
+ }
+}
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#ifndef VISU_Streamer_HeaderFile
+#define VISU_Streamer_HeaderFile
+
+#include "VISUPipeline.hxx"
+
+#include <vtkStreamer.h>
+
+class VISU_PIPELINE_EXPORT VISU_Streamer : public vtkStreamer
+{
+public:
+ vtkTypeRevisionMacro(VISU_Streamer,vtkStreamer);
+
+protected:
+ // Description:
+ // Construct object to start from position (0,0,0); integrate forward;
+ // terminal speed 0.0; vorticity computation off; integrations step length
+ // 0.2; and maximum propagation time 100.0.
+ VISU_Streamer();
+ ~VISU_Streamer();
+
+ // Integrate data
+ void Integrate(vtkDataSet *input, vtkDataSet *source);
+
+ static VTK_THREAD_RETURN_TYPE ThreadedIntegrate( void *arg );
+
+private:
+ VISU_Streamer(const VISU_Streamer&); // Not implemented.
+ void operator=(const VISU_Streamer&); // Not implemented.
+};
+
+#endif
