src/PIPELINE/VISU_StreamLine.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_StreamLine.hxx"
  21
  22 #include "vtkCellArray.h"
  23 #include "vtkDataSet.h"
  24 #include "vtkFloatArray.h"
  25 #include "vtkInformation.h"
  26 #include "vtkInformationVector.h"
  27 #include "vtkMath.h"
  28 #include "vtkObjectFactory.h"
  29 #include "vtkPointData.h"
  30 #include "vtkPolyData.h"
  31 #include "vtkPolyLine.h"
  32
  33 vtkCxxRevisionMacro(VISU_StreamLine, "$Revision$");
  34 vtkStandardNewMacro(VISU_StreamLine);
  35
  36 // Construct object with step size set to 1.0.
  37 VISU_StreamLine::VISU_StreamLine()
  38 {
  39   this->StepLength = 1.0;
  40   this->NumberOfStreamers = 0;
  41 }
  42
  43 int VISU_StreamLine::RequestData(
  44   vtkInformation *,
  45   vtkInformationVector **inputVector,
  46   vtkInformationVector *outputVector)
  47 {
  48   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  49   vtkInformation *outInfo = outputVector->GetInformationObject(0);
  50   vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
  51
  52   vtkDataSet *input = vtkDataSet::SafeDownCast(
  53     inInfo->Get(vtkDataObject::DATA_OBJECT()));
  54   vtkPolyData *output = vtkPolyData::SafeDownCast(
  55     outInfo->Get(vtkDataObject::DATA_OBJECT()));
  56   vtkDataSet *source = 0;
  57   if (sourceInfo)
  58     {
  59     source = vtkDataSet::SafeDownCast(
  60       sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
  61     }
  62
  63   vtkStreamer::StreamPoint *sPrev, *sPtr;
  64   vtkPoints *newPts;
  65   vtkFloatArray *newVectors;
  66   vtkFloatArray *newScalars=NULL;
  67   vtkCellArray *newLines;
  68   vtkIdType ptId, i, id;
  69   int j;
  70   vtkIdList *pts;
  71   double tOffset, x[3], v[3], s, r;
  72   double theta;
  73   vtkPolyLine* lineNormalGenerator = NULL;
  74   vtkFloatArray* normals = NULL;
  75   vtkFloatArray* rotation = 0;
  76
  77   this->SavePointInterval = this->StepLength;
  78   this->VISU_Streamer::Integrate(input, source);
  79   if ( this->NumberOfStreamers <= 0 ) {return 1;}
  80
  81   pts = vtkIdList::New();
  82   pts->Allocate(2500);
  83
  84   //
  85   //  Convert streamer into lines. Lines may be dashed.
  86   //
  87   newPts  = vtkPoints::New();
  88   newPts->Allocate(1000);
  89   newVectors  = vtkFloatArray::New();
  90   newVectors->SetNumberOfComponents(3);
  91   newVectors->Allocate(3000);
  92   if ( this->Vorticity )
  93     {
  94     lineNormalGenerator = vtkPolyLine::New();
  95     normals = vtkFloatArray::New();
  96     normals->SetNumberOfComponents(3);
  97     normals->Allocate(3000);
  98     rotation = vtkFloatArray::New();
  99     rotation->SetNumberOfComponents(1);
 100     rotation->Allocate(1000);
 101     rotation->SetName("Thetas");
 102     output->GetPointData()->AddArray(rotation);
 103     }
 104
 105   if ( input->GetPointData()->GetScalars() || this->SpeedScalars
 106        || this->OrientationScalars)
 107     {
 108     newScalars = vtkFloatArray::New();
 109     newScalars->Allocate(1000);
 110     }
 111   newLines = vtkCellArray::New();
 112   newLines->Allocate(newLines->EstimateSize(2*this->NumberOfStreamers,
 113                                             VTK_CELL_SIZE));
 114   //
 115   // Loop over all streamers generating points
 116   //
 117   for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
 118     {
 119     if ( this->Streamers[ptId].GetNumberOfPoints() < 2 )
 120       {
 121       continue;
 122       }
 123     sPrev = this->Streamers[ptId].GetStreamPoint(0);
 124     sPtr = this->Streamers[ptId].GetStreamPoint(1);
 125
 126     if ( this->Streamers[ptId].GetNumberOfPoints() == 2 && sPtr->cellId >= 0 )
 127       {
 128       continue;
 129       }
 130
 131     tOffset = sPrev->t;
 132
 133     for ( i=1;
 134     i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
 135     i++, sPrev=sPtr, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
 136       {
 137       //
 138       // Create points for line
 139       //
 140       while ( tOffset >= sPrev->t && tOffset < sPtr->t )
 141         {
 142         r = (tOffset - sPrev->t) / (sPtr->t - sPrev->t);
 143
 144         for (j=0; j<3; j++)
 145           {
 146           x[j] = sPrev->x[j] + r * (sPtr->x[j] - sPrev->x[j]);
 147           v[j] = sPrev->v[j] + r * (sPtr->v[j] - sPrev->v[j]);
 148           }
 149
 150         // add point to line
 151         id = newPts->InsertNextPoint(x);
 152         pts->InsertNextId(id);
 153         newVectors->InsertTuple(id,v);
 154
 155         if ( newScalars )
 156           {
 157           s = sPrev->s + r * (sPtr->s - sPrev->s);
 158           newScalars->InsertTuple(id,&s);
 159           }
 160
 161         if ( this->Vorticity )
 162           {
 163           // Store the rotation values. Used after all the streamlines
 164           // are generated.
 165           theta = sPrev->theta + r * (sPtr->theta - sPrev->theta);
 166           rotation->InsertTuple(id, &theta);
 167           }
 168
 169         tOffset += this->StepLength;
 170
 171         } // while
 172       } //for this streamer
 173
 174     if ( pts->GetNumberOfIds() > 1 )
 175       {
 176       newLines->InsertNextCell(pts);
 177       pts->Reset();
 178       }
 179     } //for all streamers
 180
 181   vtkDebugMacro(<<"Created " << newPts->GetNumberOfPoints() << " points, "
 182                << newLines->GetNumberOfCells() << " lines");
 183
 184   if (this->Vorticity)
 185     {
 186     // Rotate the normal vectors with stream vorticity
 187     vtkIdType nPts=0;
 188     vtkIdType *linePts=0;
 189     double normal[3], local1[3], local2[3], length, costheta, sintheta;
 190
 191     lineNormalGenerator->GenerateSlidingNormals(newPts,newLines,normals);
 192
 193     //  Loop over all lines, from the above code we are know that each line
 194     //  will have at least two points and that no points will be shared
 195     //  between lines. It is important to loop over the points used by the
 196     //  lines because newPts may actually contain points that are not used by
 197     //  any lines. The normals are only calculated for points that are used
 198     //  in lines so referencing normals for all points can lead to UMRs
 199     for (newLines->InitTraversal(); newLines->GetNextCell(nPts,linePts); )
 200       {
 201       for(i=0; i<nPts; i++)
 202         {
 203         normals->GetTuple(linePts[i], normal);
 204         newVectors->GetTuple(linePts[i], v);
 205         // obtain two unit orthogonal vectors on the plane perpendicular to
 206         // the streamline
 207         for(j=0; j<3; j++)
 208           {
 209           local1[j] = normal[j];
 210           }
 211         length = vtkMath::Normalize(local1);
 212         vtkMath::Cross(local1, v, local2);
 213         vtkMath::Normalize(local2);
 214         // Rotate the normal with theta
 215         rotation->GetTuple(linePts[i], &theta);
 216         costheta = cos(theta);
 217         sintheta = sin(theta);
 218         for(j=0; j<3; j++)
 219           {
 220           normal[j] = length* (costheta*local1[j] + sintheta*local2[j]);
 221           }
 222         normals->SetTuple(linePts[i], normal);
 223         }
 224       }
 225     output->GetPointData()->SetNormals(normals);
 226     normals->Delete();
 227     lineNormalGenerator->Delete();
 228     rotation->Delete();
 229     }
 230
 231   output->SetPoints(newPts);
 232   newPts->Delete();
 233
 234   output->GetPointData()->SetVectors(newVectors);
 235   newVectors->Delete();
 236
 237   if ( newScalars )
 238     {
 239     int idx = output->GetPointData()->AddArray(newScalars);
 240     output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
 241     newScalars->Delete();
 242     }
 243
 244   pts->Delete();
 245   output->SetLines(newLines);
 246   newLines->Delete();
 247
 248   // Delete the streamers since they are no longer needed
 249   delete[] this->Streamers;
 250   this->Streamers = 0;
 251   this->NumberOfStreamers = 0;
 252
 253   output->Squeeze();
 254
 255   return 1;
 256 }
 257
 258 void VISU_StreamLine::PrintSelf(ostream& os, vtkIndent indent)
 259 {
 260   this->Superclass::PrintSelf(os,indent);
 261
 262   os << indent << "Step Length: " << this->StepLength << "\n";
 263
 264 }