src/Plugins/DifferenceTimesteps/vtkDifferenceTimestepsFilter.cxx

   1 // Copyright (C) 2014  CEA/DEN, EDF R&D
   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, or (at your option) any later version.
   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 // Author : Maxim Glibin
  20
  21 #include "vtkDifferenceTimestepsFilter.h"
  22
  23 #include <vtkCellData.h>
  24 #include <vtkDataArray.h>
  25 #include <vtkDataObjectTreeIterator.h>
  26 #include <vtkDataSet.h>
  27 #include <vtkDoubleArray.h>
  28 #include <vtkInformation.h>
  29 #include <vtkInformationVector.h>
  30 #include <vtkMultiBlockDataSet.h>
  31 #include <vtkObjectFactory.h>
  32 #include <vtkPointData.h>
  33 #include <vtkStreamingDemandDrivenPipeline.h>
  34 #include <vtkStringArray.h>
  35 #include <vtkUnstructuredGrid.h>
  36
  37 // Temporal difference of data array
  38 vtkDataArray* DataTempDiffArray(vtkDataArray* theDataArray,
  39                                 vtkIdType     theNumComp,
  40                                 vtkIdType     theNumTuple,
  41                                 const char*   thePrefix)
  42 {
  43   // Create the new array
  44   vtkAbstractArray *anAbstractArray = theDataArray->CreateArray(theDataArray->GetDataType());
  45   vtkDataArray *anOutput = vtkDataArray::SafeDownCast(anAbstractArray);
  46
  47   // Initialize and appoint a new name
  48   anOutput->SetNumberOfComponents(theNumComp);
  49   anOutput->SetNumberOfTuples(theNumTuple);
  50   vtkstd::string aNewName = vtkstd::string(thePrefix) + theDataArray->GetName();
  51   anOutput->SetName(aNewName.c_str());
  52
  53   return anOutput;
  54 }
  55
  56 // Templated difference function
  57 template <class T>
  58 void vtkTemporalDataDifference(
  59   vtkDifferenceTimestepsFilter* theDTF,
  60   vtkDataArray*                 theOutput,
  61   vtkDataArray**                theArrays,
  62   vtkIdType                     theNumComp,
  63   T *)
  64 {
  65   T* anOutputData = static_cast<T*>(theOutput->GetVoidPointer(0));
  66   T* anInputData0 = static_cast<T*>(theArrays[0]->GetVoidPointer(0));
  67   T* anInputData1 = static_cast<T*>(theArrays[1]->GetVoidPointer(0));
  68
  69   vtkIdType N = theArrays[0]->GetNumberOfTuples();
  70   for (vtkIdType t = 0; t < N; ++t)
  71   {
  72     T* x0 = &anInputData0[t*theNumComp];
  73     T* x1 = &anInputData1[t*theNumComp];
  74     for (int c = 0; c < theNumComp; ++c)
  75     {
  76       // Compute the difference
  77       *anOutputData++ = static_cast<T>(x1[c]-x0[c]);
  78     }
  79   }
  80   theOutput->SetNumberOfTuples(N);
  81 }
  82
  83 vtkStandardNewMacro(vtkDifferenceTimestepsFilter);
  84
  85 //--------------------------------------------------------------------------------------------------
  86 vtkDifferenceTimestepsFilter::vtkDifferenceTimestepsFilter()
  87 {
  88   this->NumberTimeSteps = 0;
  89   this->RangeIndicesTimeSteps[0] = 0;
  90   this->RangeIndicesTimeSteps[1] = 0;
  91   this->FirstTimeStepIndex = 0.0;
  92   this->SecondTimeStepIndex = 0.0;
  93   this->TimeStepValues.clear();
  94   this->ArrayNamePrefix = NULL;
  95
  96   this->SetNumberOfInputPorts(1);
  97   this->SetNumberOfOutputPorts(1);
  98
  99   // Set the input data array that the algorithm will process
 100   this->SetInputArrayToProcess(
 101     0,
 102     0,
 103     0,
 104     vtkDataObject::FIELD_ASSOCIATION_POINTS,
 105     vtkDataSetAttributes::SCALARS);
 106 }
 107
 108 //--------------------------------------------------------------------------------------------------
 109 vtkDifferenceTimestepsFilter::~vtkDifferenceTimestepsFilter()
 110 {
 111   this->TimeStepValues.clear();
 112   this->SetArrayNamePrefix(NULL);
 113 }
 114
 115 //--------------------------------------------------------------------------------------------------
 116 void vtkDifferenceTimestepsFilter::PrintSelf(ostream& theOS, vtkIndent theIndent)
 117 {
 118   this->Superclass::PrintSelf(theOS, theIndent);
 119   theOS << theIndent << "Number of time steps : " << this->NumberTimeSteps << endl;
 120   theOS << theIndent << "First time step : " << this->FirstTimeStepIndex << endl;
 121   theOS << theIndent << "Second time step : " << this->SecondTimeStepIndex << endl;
 122   theOS << theIndent << "Field association : "
 123         << vtkDataObject::GetAssociationTypeAsString(this->GetInputFieldAssociation()) << endl;
 124 }
 125
 126 //--------------------------------------------------------------------------------------------------
 127 int vtkDifferenceTimestepsFilter::FillInputPortInformation(
 128   int thePort, vtkInformation* theInfo)
 129 {
 130   if (thePort == 0)
 131     theInfo->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
 132
 133   return 1;
 134 }
 135
 136 //--------------------------------------------------------------------------------------------------
 137 int vtkDifferenceTimestepsFilter::FillOutputPortInformation(
 138   int vtkNotUsed(thePort), vtkInformation* theInfo)
 139 {
 140   theInfo->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkDataObject");
 141   return 1;
 142 }
 143
 144 //--------------------------------------------------------------------------------------------------
 145 int vtkDifferenceTimestepsFilter::RequestDataObject(
 146   vtkInformation*        vtkNotUsed(theRequest),
 147   vtkInformationVector** theInputVector,
 148   vtkInformationVector*  theOutputVector)
 149 {
 150   if (this->GetNumberOfInputPorts() == 0 || this->GetNumberOfOutputPorts() == 0)
 151     return 1;
 152
 153   vtkInformation* anInputInfo = theInputVector[0]->GetInformationObject(0);
 154   if (anInputInfo == NULL)
 155   {
 156     vtkErrorMacro(<< "Input information vector is missed.");
 157     return 0;
 158   }
 159
 160   vtkDataObject *anInputObj = anInputInfo->Get(vtkDataObject::DATA_OBJECT());
 161   if (anInputObj != NULL)
 162   {
 163     // For each output
 164     for (int i = 0; i < this->GetNumberOfOutputPorts(); ++i)
 165     {
 166       vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(i);
 167       vtkDataObject *anOutputObj = anOutputInfo->Get(vtkDataObject::DATA_OBJECT());
 168       if (!anOutputObj || !anOutputObj->IsA(anInputObj->GetClassName()))
 169       {
 170         vtkDataObject* aNewOutput = anInputObj->NewInstance();
 171         anOutputInfo->Set(vtkDataObject::DATA_OBJECT(), aNewOutput);
 172         aNewOutput->Delete();
 173       }
 174     }
 175     return 1;
 176   }
 177   return 0;
 178 }
 179
 180 //--------------------------------------------------------------------------------------------------
 181 int vtkDifferenceTimestepsFilter::RequestInformation(
 182   vtkInformation*        vtkNotUsed(theRequest),
 183   vtkInformationVector** theInputVector,
 184   vtkInformationVector*  theOutputVector)
 185 {
 186   // Get input and output information objects
 187   vtkInformation *anInInfo = theInputVector[0]->GetInformationObject(0);
 188   vtkInformation *anOutInfo = theOutputVector->GetInformationObject(0);
 189
 190   // Check for presence more than one time step
 191   if (anInInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
 192   {
 193     // Find time on input
 194     this->NumberTimeSteps = anInInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
 195     if (this->NumberTimeSteps < 2)
 196     {
 197       vtkErrorMacro(<< "Not enough numbers of time steps: " << this->NumberTimeSteps);
 198       return 0;
 199     }
 200     // Get time step values
 201     this->TimeStepValues.resize(this->NumberTimeSteps);
 202     anInInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &this->TimeStepValues[0]);
 203     if (this->TimeStepValues.size() == 0)
 204     {
 205       vtkErrorMacro(<<"Array of time steps is empty.");
 206       return 0;
 207     }
 208   }
 209   else
 210   {
 211     vtkErrorMacro(<< "No time steps in input data.");
 212     return 0;
 213   }
 214
 215   // Update range of indices of the time steps
 216   this->RangeIndicesTimeSteps[0] = 0;
 217   this->RangeIndicesTimeSteps[1] = this->NumberTimeSteps - 1;
 218
 219   // The output data of this filter has no time associated with it.
 220   // It is the result of computation difference between two time steps.
 221   // Unset the time steps
 222   if (anOutInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
 223     anOutInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
 224
 225   // Unset the time range
 226   if(anOutInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_RANGE()))
 227     anOutInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
 228
 229   return 1;
 230 }
 231
 232 //--------------------------------------------------------------------------------------------------
 233 int vtkDifferenceTimestepsFilter::RequestUpdateExtent(vtkInformation*        theRequest,
 234                                                       vtkInformationVector** theInputVector,
 235                                                       vtkInformationVector*  theOutputVector)
 236 {
 237   // Get the information objects
 238   vtkInformation* anInputInfo  = theInputVector[0]->GetInformationObject(0);
 239   vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(0);
 240
 241   // Indices must not go beyond the range of indices of the time steps
 242   if (this->FirstTimeStepIndex >= this->NumberTimeSteps || this->FirstTimeStepIndex < 0)
 243   {
 244     vtkErrorMacro(<< "Specified index of the first time step ["
 245                   << this->FirstTimeStepIndex
 246                   << "] is outside the range of indices.");
 247     return 0;
 248   }
 249   if (this->SecondTimeStepIndex >= this->NumberTimeSteps || this->SecondTimeStepIndex < 0)
 250   {
 251     vtkErrorMacro(<< "Specified index of the second time step ["
 252                   << this->SecondTimeStepIndex
 253                   << "] is outside the range of indices.");
 254     return 0;
 255   }
 256
 257   // Warn if the selected time steps are equal
 258   if (this->FirstTimeStepIndex == this->SecondTimeStepIndex)
 259   {
 260     vtkWarningMacro(<< "First and second indices ["
 261                     << this->FirstTimeStepIndex
 262                     << " = " << this->SecondTimeStepIndex
 263                     << "] are the same.");
 264   }
 265
 266   // Find the required input time steps and request them
 267   if (anOutputInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
 268   {
 269     // Get the available input times
 270     double *anInputTimes = anInputInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
 271     if (anInputTimes != NULL)
 272     {
 273       // Compute the requested times
 274       double anInputUpdateTimes[2];
 275       int aNumberInputUpdateTimes(0);
 276
 277       // For each the requested time mark the required input times
 278       anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->FirstTimeStepIndex];
 279       anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->SecondTimeStepIndex];
 280
 281       // Make the multiple time requests upstream and use set of time-stamped data
 282       // objects are stored in time order in a vtkMultiBlockDataSet object
 283       anInputInfo->Set(vtkMultiTimeStepAlgorithm::UPDATE_TIME_STEPS(),
 284                        anInputUpdateTimes, aNumberInputUpdateTimes);
 285     }
 286   }
 287   return 1;
 288 }
 289
 290 //--------------------------------------------------------------------------------------------------
 291 int vtkDifferenceTimestepsFilter::RequestData(vtkInformation*        vtkNotUsed(theRequest),
 292                                               vtkInformationVector** theInputVector,
 293                                               vtkInformationVector*  theOutputVector)
 294 {
 295   // Get the information objects
 296   vtkInformation *anInputInfo  = theInputVector[0]->GetInformationObject(0);
 297   vtkInformation *anOutputInfo = theOutputVector->GetInformationObject(0);
 298
 299   vtkDataObject *anOutputDataObj = NULL;
 300
 301   vtkMultiBlockDataSet *anInputData =
 302     vtkMultiBlockDataSet::SafeDownCast(anInputInfo->Get(vtkDataObject::DATA_OBJECT()));
 303
 304   int aNumberTimeSteps = anInputData->GetNumberOfBlocks();
 305   if (aNumberTimeSteps == 2)
 306   {
 307     // Get data objects
 308     vtkDataObject* aData0 = anInputData->GetBlock(0);
 309     vtkDataObject* aData1 = anInputData->GetBlock(1);
 310     if (aData0 == NULL && aData1 == NULL)
 311     {
 312       vtkErrorMacro(<< "Null data set.");
 313       return 0;
 314     }
 315
 316     // Compute difference between two objects
 317     anOutputDataObj = this->DifferenceDataObject(aData0, aData1);
 318     anOutputInfo->Set(vtkDataObject::DATA_OBJECT(), anOutputDataObj);
 319     if (anOutputDataObj != NULL)
 320       anOutputDataObj->Delete();
 321   }
 322   else
 323   {
 324     vtkErrorMacro(<< "The amount of time blocks is not correct: " << aNumberTimeSteps);
 325     return 0;
 326   }
 327
 328   return 1;
 329 }
 330
 331 //--------------------------------------------------------------------------------------------------
 332 vtkDataObject* vtkDifferenceTimestepsFilter::DifferenceDataObject(vtkDataObject* theInput1,
 333                                                                   vtkDataObject* theInput2)
 334 {
 335   // Determine the input object type
 336   if (theInput1->IsA("vtkDataSet"))
 337   {
 338     vtkDataSet *anInDataSet1 = vtkDataSet::SafeDownCast(theInput1);
 339     vtkDataSet *anInDataSet2 = vtkDataSet::SafeDownCast(theInput2);
 340     return this->DifferenceDataSet(anInDataSet1, anInDataSet2);
 341   }
 342   else if (theInput1->IsA("vtkCompositeDataSet"))
 343   {
 344     // It is essential that aMGDataSet[0] an aMGDataSet[1] has the same structure.
 345     vtkCompositeDataSet* aMGDataSet[2];
 346     aMGDataSet[0] = vtkCompositeDataSet::SafeDownCast(theInput1);
 347     aMGDataSet[1] = vtkCompositeDataSet::SafeDownCast(theInput2);
 348
 349     vtkCompositeDataSet *anOutput = aMGDataSet[0]->NewInstance();
 350     anOutput->CopyStructure(aMGDataSet[0]);
 351
 352     vtkSmartPointer<vtkCompositeDataIterator> anIter;
 353     anIter.TakeReference(aMGDataSet[0]->NewIterator());
 354     for (anIter->InitTraversal(); !anIter->IsDoneWithTraversal(); anIter->GoToNextItem())
 355     {
 356       vtkDataObject* aDataObj1 = anIter->GetCurrentDataObject();
 357       vtkDataObject* aDataObj2 = aMGDataSet[1]->GetDataSet(anIter);
 358       if (aDataObj1 == NULL || aDataObj2 == NULL)
 359       {
 360         vtkWarningMacro("The composite datasets were not identical in structure.");
 361         continue;
 362       }
 363
 364       vtkDataObject *aResultDObj = this->DifferenceDataObject(aDataObj1, aDataObj2);
 365       if (aResultDObj != NULL)
 366       {
 367         anOutput->SetDataSet(anIter, aResultDObj);
 368         aResultDObj->Delete();
 369       }
 370       else
 371       {
 372         vtkErrorMacro(<< "Unexpected error during computation of the difference.");
 373         return NULL;
 374       }
 375     }
 376     return anOutput;
 377   }
 378   else
 379   {
 380     vtkErrorMacro("We cannot yet compute difference of this type of dataset.");
 381     return NULL;
 382   }
 383 }
 384
 385 //--------------------------------------------------------------------------------------------------
 386 vtkDataSet* vtkDifferenceTimestepsFilter::DifferenceDataSet(vtkDataSet* theInput1,
 387                                                             vtkDataSet* theInput2)
 388 {
 389   vtkDataSet *anInput[2];
 390   anInput[0] = theInput1;
 391   anInput[1] = theInput2;
 392
 393   // Copy input structure into output
 394   vtkDataSet *anOutput = anInput[0]->NewInstance();
 395   anOutput->CopyStructure(anInput[0]);
 396
 397   vtkstd::vector<vtkDataArray*> anArrays;
 398   vtkDataArray *anOutputArray;
 399
 400   // Compute the difference of the the specified point or cell data array
 401   vtkDataArray* aDataArray0 = this->GetInputArrayToProcess(0, anInput[0]);
 402   vtkDataArray* aDataArray1 = this->GetInputArrayToProcess(0, anInput[1]);
 403   if (aDataArray0 == NULL || aDataArray1 == NULL)
 404   {
 405     vtkErrorMacro(<< "Input array to process is empty.");
 406     return NULL;
 407   }
 408   anArrays.push_back(aDataArray0);
 409   anArrays.push_back(aDataArray1);
 410
 411   if (anArrays.size() > 1)
 412   {
 413     if (!this->VerifyArrays(&anArrays[0], 2))
 414     {
 415       vtkErrorMacro(<< "Verification of data arrays has failed.");
 416       return NULL;
 417     }
 418
 419     anOutputArray = this->DifferenceDataArray(&anArrays[0], anArrays[0]->GetNumberOfTuples());
 420     // Determine a field association
 421     if (this->GetInputFieldAssociation() == vtkDataObject::FIELD_ASSOCIATION_POINTS)
 422     {
 423       // For point data
 424       anOutput->GetPointData()->AddArray(anOutputArray);
 425     }
 426     else if (this->GetInputFieldAssociation() == vtkDataObject::FIELD_ASSOCIATION_CELLS)
 427     {
 428       // For cell data
 429       anOutput->GetCellData()->AddArray(anOutputArray);
 430     }
 431     else
 432     {
 433       vtkErrorMacro(<< "Solution is not implemeted yet.");
 434       return NULL;
 435     }
 436     anOutputArray->Delete();
 437     anArrays.clear();
 438   }
 439
 440   return anOutput;
 441 }
 442
 443 //--------------------------------------------------------------------------------------------------
 444 vtkDataArray* vtkDifferenceTimestepsFilter::DifferenceDataArray(vtkDataArray** theArrays,
 445                                                                 vtkIdType      theNumTuple)
 446 {
 447   // Create the output array based on the number of tuple and components
 448   // with a new name containing the specified prefix
 449   int aNumComp = theArrays[0]->GetNumberOfComponents();
 450   vtkDataArray *anOutput =
 451     DataTempDiffArray(theArrays[0], aNumComp, theNumTuple, this->ArrayNamePrefix);
 452
 453   // Now do the computation of the difference
 454   switch (theArrays[0]->GetDataType())
 455   {
 456     vtkTemplateMacro(
 457       vtkTemporalDataDifference(this, anOutput, theArrays, aNumComp, static_cast<VTK_TT *>(0)));
 458   default:
 459     vtkErrorMacro(<< "Execute: unknown scalar type.");
 460     return NULL;
 461   }
 462
 463   return anOutput;
 464 }
 465
 466 //--------------------------------------------------------------------------------------------------
 467 int vtkDifferenceTimestepsFilter::GetInputFieldAssociation()
 468 {
 469   vtkInformationVector *anInputArrayVec = this->GetInformation()->Get(INPUT_ARRAYS_TO_PROCESS());
 470   vtkInformation *anInputArrayInfo = anInputArrayVec->GetInformationObject(0);
 471   return anInputArrayInfo->Get(vtkDataObject::FIELD_ASSOCIATION());
 472 }
 473
 474 //--------------------------------------------------------------------------------------------------
 475 bool vtkDifferenceTimestepsFilter::VerifyArrays(vtkDataArray **theArrays, int theNumArrays)
 476 {
 477   // Get all required data to compare with other
 478   const char* anArrayName  = theArrays[0]->GetName();
 479   vtkIdType aNumTuples     = theArrays[0]->GetNumberOfTuples();
 480   vtkIdType aNumComponents = theArrays[0]->GetNumberOfComponents();
 481
 482   for (int i = 1; i < theNumArrays; ++i)
 483   {
 484     if (strcmp(theArrays[i]->GetName(), anArrayName) != 0)
 485     {
 486       vtkWarningMacro(<< "Computation of difference aborted for dataset because "
 487                       << "the array name in each time step are different.")
 488       return false;
 489     }
 490
 491     if (theArrays[i]->GetNumberOfTuples() != aNumTuples)
 492     {
 493       vtkWarningMacro(<< "Computation of difference aborted for dataset because "
 494                       << "the number of tuples in each time step are different.")
 495       return false;
 496     }
 497
 498     if (theArrays[i]->GetNumberOfComponents() != aNumComponents)
 499     {
 500       vtkWarningMacro(<< "Computation of difference aborted for dataset because "
 501                       << "the number of components in each time step are different.")
 502       return false;
 503     }
 504   }
 505
 506   return true;
 507 }