From: ouv Date: Tue, 9 Mar 2010 10:11:17 +0000 (+0000) Subject: Issue 0020724: EDF 1237 VISU : Crash Salome with Stream Lines X-Git-Tag: V5_1_main_20100310 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=6fff75f453e0fbc8ba9d0bd66e845fb922d5d17a;p=modules%2Fvisu.git Issue 0020724: EDF 1237 VISU : Crash Salome with Stream Lines --- diff --git a/src/PIPELINE/Makefile.am b/src/PIPELINE/Makefile.am index 1ef1a09e..a7558e34 100644 --- a/src/PIPELINE/Makefile.am +++ b/src/PIPELINE/Makefile.am @@ -46,6 +46,8 @@ salomeinclude_HEADERS= \ VISU_IsoSurfacesPL.hxx \ VISU_DeformedShapePL.hxx \ VISU_VectorsPL.hxx \ + VISU_Streamer.hxx \ + VISU_StreamLine.hxx \ VISU_StreamLinesPL.hxx \ VISU_LookupTable.hxx \ VISU_ScalarBarActor.hxx \ @@ -92,6 +94,8 @@ dist_libVisuPipeLine_la_SOURCES= \ VISU_IsoSurfacesPL.cxx \ VISU_DeformedShapePL.cxx \ VISU_VectorsPL.cxx \ + VISU_Streamer.cxx \ + VISU_StreamLine.cxx \ VISU_StreamLinesPL.cxx \ VISU_LookupTable.cxx \ VISU_ScalarBarActor.cxx \ diff --git a/src/PIPELINE/VISU_PipeLine.cxx b/src/PIPELINE/VISU_PipeLine.cxx index 4ffd2030..89bf6b40 100644 --- a/src/PIPELINE/VISU_PipeLine.cxx +++ b/src/PIPELINE/VISU_PipeLine.cxx @@ -444,8 +444,8 @@ VISU_PipeLine //---------------------------------------------------------------------------- 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; diff --git a/src/PIPELINE/VISU_PipeLine.hxx b/src/PIPELINE/VISU_PipeLine.hxx index ebb73c7e..959142dd 100644 --- a/src/PIPELINE/VISU_PipeLine.hxx +++ b/src/PIPELINE/VISU_PipeLine.hxx @@ -185,8 +185,8 @@ public: static size_t - GetAvailableMemory(size_t theSize, - size_t theMinSize = 1024*1024); + GetAvailableMemory(double theSize, + double theMinSize = 1024*1024); protected: //---------------------------------------------------------------------------- diff --git a/src/PIPELINE/VISU_StreamLine.cxx b/src/PIPELINE/VISU_StreamLine.cxx new file mode 100644 index 00000000..9c0979fd --- /dev/null +++ b/src/PIPELINE/VISU_StreamLine.cxx @@ -0,0 +1,266 @@ +// 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; iGetTuple(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"; + +} diff --git a/src/PIPELINE/VISU_StreamLine.hxx b/src/PIPELINE/VISU_StreamLine.hxx new file mode 100644 index 00000000..150094a3 --- /dev/null +++ b/src/PIPELINE/VISU_StreamLine.hxx @@ -0,0 +1,60 @@ +// 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 diff --git a/src/PIPELINE/VISU_StreamLinesPL.cxx b/src/PIPELINE/VISU_StreamLinesPL.cxx index 1a832266..40d0dcd8 100644 --- a/src/PIPELINE/VISU_StreamLinesPL.cxx +++ b/src/PIPELINE/VISU_StreamLinesPL.cxx @@ -31,6 +31,7 @@ //#include "VISU_UsedPointsFilter.hxx" #include "VISU_MaskPointsFilter.hxx" #include "VISU_PipeLineUtils.hxx" +#include "VISU_StreamLine.hxx" #include "VTKViewer_GeometryFilter.h" @@ -38,7 +39,6 @@ #include #include -#include #ifdef _DEBUG_ static int MYDEBUG = 0; @@ -63,7 +63,7 @@ VISU_StreamLinesPL SetIsShrinkable(false); SetIsFeatureEdgesAllowed(false); - myStream = vtkStreamLine::New(); + myStream = VISU_StreamLine::New(); myCenters = vtkCellCenters::New(); myGeomFilter = VTKViewer_GeometryFilter::New(); myPointsFilter = VISU_MaskPointsFilter::New(); @@ -129,7 +129,7 @@ VISU_StreamLinesPL //---------------------------------------------------------------------------- -size_t +vtkFloatingPointType VISU_StreamLinesPL ::GetNecasseryMemorySize(vtkIdType theNbOfPoints, vtkFloatingPointType theStepLength, @@ -153,7 +153,7 @@ VISU_StreamLinesPL 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; } @@ -168,7 +168,7 @@ VISU_StreamLinesPL { 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); @@ -389,10 +389,11 @@ VISU_StreamLinesPL 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))); diff --git a/src/PIPELINE/VISU_StreamLinesPL.hxx b/src/PIPELINE/VISU_StreamLinesPL.hxx index d2ec3a86..fd72d0f3 100644 --- a/src/PIPELINE/VISU_StreamLinesPL.hxx +++ b/src/PIPELINE/VISU_StreamLinesPL.hxx @@ -29,13 +29,15 @@ #include "VISUPipeline.hxx" #include "VISU_DeformedShapePL.hxx" -#include + +#include class vtkDataSet; class vtkPointSet; class vtkCellCenters; class VTKViewer_GeometryFilter; class VISU_MaskPointsFilter; +class VISU_StreamLine; //---------------------------------------------------------------------------- @@ -207,7 +209,7 @@ protected: bool theIsCopyInput); static - size_t + vtkFloatingPointType GetNecasseryMemorySize(vtkIdType theNbOfPoints, vtkFloatingPointType theStepLength, vtkFloatingPointType thePropogationTime, @@ -236,7 +238,7 @@ protected: CorrectStepLength(vtkFloatingPointType theStep, vtkDataSet* theDataSet); - vtkStreamLine* myStream; + VISU_StreamLine* myStream; vtkPointSet* mySource; vtkCellCenters* myCenters; VTKViewer_GeometryFilter *myGeomFilter; diff --git a/src/PIPELINE/VISU_Streamer.cxx b/src/PIPELINE/VISU_Streamer.cxx new file mode 100644 index 00000000..6f8cbc89 --- /dev/null +++ b/src/PIPELINE/VISU_Streamer.cxx @@ -0,0 +1,571 @@ +// 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(); + } +} diff --git a/src/PIPELINE/VISU_Streamer.hxx b/src/PIPELINE/VISU_Streamer.hxx new file mode 100644 index 00000000..987a8973 --- /dev/null +++ b/src/PIPELINE/VISU_Streamer.hxx @@ -0,0 +1,52 @@ +// 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 + +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