src/Plugins/DifferenceTimesteps/vtkDifferenceTimestepsFilter.cxx

   1 // Copyright (C) 2014-2019  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   std::string aNewName = std::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   /*
 220    * RNV: Temporary commented:
 221    *      This piece of the code removes all time steps from the output object,
 222    *      but this leads to the strange side effect in the ParaView: time steps also disappears
 223    *      from the animation scene of the input (parent) object of this filter.
 224    *      Seems it is a bug of the ParaView, to be investigated ...
 225    *
 226   // The output data of this filter has no time associated with it.
 227   // It is the result of computation difference between two time steps.
 228   // Unset the time steps
 229   if (anOutInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
 230     anOutInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
 231
 232   // Unset the time range
 233   if(anOutInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_RANGE()))
 234     anOutInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
 235   */
 236   return 1;
 237 }
 238
 239 //--------------------------------------------------------------------------------------------------
 240 int vtkDifferenceTimestepsFilter::RequestUpdateExtent(vtkInformation*        theRequest,
 241                                                       vtkInformationVector** theInputVector,
 242                                                       vtkInformationVector*  theOutputVector)
 243 {
 244   // Get the information objects
 245   vtkInformation* anInputInfo  = theInputVector[0]->GetInformationObject(0);
 246   vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(0);
 247
 248   // Indices must not go beyond the range of indices of the time steps
 249   if (this->FirstTimeStepIndex >= this->NumberTimeSteps || this->FirstTimeStepIndex < 0)
 250   {
 251     vtkErrorMacro(<< "Specified index of the first time step ["
 252                   << this->FirstTimeStepIndex
 253                   << "] is outside the range of indices.");
 254     return 0;
 255   }
 256   if (this->SecondTimeStepIndex >= this->NumberTimeSteps || this->SecondTimeStepIndex < 0)
 257   {
 258     vtkErrorMacro(<< "Specified index of the second time step ["
 259                   << this->SecondTimeStepIndex
 260                   << "] is outside the range of indices.");
 261     return 0;
 262   }
 263
 264   // Warn if the selected time steps are equal
 265   if (this->FirstTimeStepIndex == this->SecondTimeStepIndex)
 266   {
 267     vtkWarningMacro(<< "First and second indices ["
 268                     << this->FirstTimeStepIndex
 269                     << " = " << this->SecondTimeStepIndex
 270                     << "] are the same.");
 271   }
 272
 273   // Find the required input time steps and request them
 274   if (anOutputInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
 275   {
 276     // Get the available input times
 277     double *anInputTimes = anInputInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
 278     if (anInputTimes != NULL)
 279     {
 280       // Compute the requested times
 281       double anInputUpdateTimes[2];
 282       int aNumberInputUpdateTimes(0);
 283
 284       // For each the requested time mark the required input times
 285       anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->FirstTimeStepIndex];
 286       anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->SecondTimeStepIndex];
 287
 288       // Make the multiple time requests upstream and use set of time-stamped data
 289       // objects are stored in time order in a vtkMultiBlockDataSet object
 290       anInputInfo->Set(vtkMultiTimeStepAlgorithm::UPDATE_TIME_STEPS(),
 291                        anInputUpdateTimes, aNumberInputUpdateTimes);
 292     }
 293   }
 294   return 1;
 295 }
 296
 297 //--------------------------------------------------------------------------------------------------
 298 int vtkDifferenceTimestepsFilter::RequestData(vtkInformation*        vtkNotUsed(theRequest),
 299                                               vtkInformationVector** theInputVector,
 300                                               vtkInformationVector*  theOutputVector)
 301 {
 302   // Get the information objects
 303   vtkInformation *anInputInfo  = theInputVector[0]->GetInformationObject(0);
 304   vtkInformation *anOutputInfo = theOutputVector->GetInformationObject(0);
 305
 306   vtkDataObject *anOutputDataObj = NULL;
 307
 308   vtkMultiBlockDataSet *anInputData =
 309     vtkMultiBlockDataSet::SafeDownCast(anInputInfo->Get(vtkDataObject::DATA_OBJECT()));
 310
 311   int aNumberTimeSteps = anInputData->GetNumberOfBlocks();
 312   if (aNumberTimeSteps == 2)
 313   {
 314     // Get data objects
 315     vtkDataObject* aData0 = anInputData->GetBlock(0);
 316     vtkDataObject* aData1 = anInputData->GetBlock(1);
 317     if (aData0 == NULL && aData1 == NULL)
 318     {
 319       vtkErrorMacro(<< "Null data set.");
 320       return 0;
 321     }
 322
 323     // Compute difference between two objects
 324     anOutputDataObj = this->DifferenceDataObject(aData0, aData1);
 325     anOutputInfo->Set(vtkDataObject::DATA_OBJECT(), anOutputDataObj);
 326     if (anOutputDataObj != NULL)
 327       anOutputDataObj->Delete();
 328   }
 329   else
 330   {
 331     vtkErrorMacro(<< "The amount of time blocks is not correct: " << aNumberTimeSteps);
 332     return 0;
 333   }
 334
 335   return 1;
 336 }
 337
 338 //--------------------------------------------------------------------------------------------------
 339 vtkDataObject* vtkDifferenceTimestepsFilter::DifferenceDataObject(vtkDataObject* theInput1,
 340                                                                   vtkDataObject* theInput2)
 341 {
 342   // Determine the input object type
 343   if (theInput1->IsA("vtkDataSet"))
 344   {
 345     vtkDataSet *anInDataSet1 = vtkDataSet::SafeDownCast(theInput1);
 346     vtkDataSet *anInDataSet2 = vtkDataSet::SafeDownCast(theInput2);
 347     return this->DifferenceDataSet(anInDataSet1, anInDataSet2);
 348   }
 349   else if (theInput1->IsA("vtkCompositeDataSet"))
 350   {
 351     // It is essential that aMGDataSet[0] an aMGDataSet[1] has the same structure.
 352     vtkCompositeDataSet* aMGDataSet[2];
 353     aMGDataSet[0] = vtkCompositeDataSet::SafeDownCast(theInput1);
 354     aMGDataSet[1] = vtkCompositeDataSet::SafeDownCast(theInput2);
 355
 356     vtkCompositeDataSet *anOutput = aMGDataSet[0]->NewInstance();
 357     anOutput->CopyStructure(aMGDataSet[0]);
 358
 359     vtkSmartPointer<vtkCompositeDataIterator> anIter;
 360     anIter.TakeReference(aMGDataSet[0]->NewIterator());
 361     for (anIter->InitTraversal(); !anIter->IsDoneWithTraversal(); anIter->GoToNextItem())
 362     {
 363       vtkDataObject* aDataObj1 = anIter->GetCurrentDataObject();
 364       vtkDataObject* aDataObj2 = aMGDataSet[1]->GetDataSet(anIter);
 365       if (aDataObj1 == NULL || aDataObj2 == NULL)
 366       {
 367         vtkWarningMacro("The composite datasets were not identical in structure.");
 368         continue;
 369       }
 370
 371       vtkDataObject *aResultDObj = this->DifferenceDataObject(aDataObj1, aDataObj2);
 372       if (aResultDObj != NULL)
 373       {
 374         anOutput->SetDataSet(anIter, aResultDObj);
 375         aResultDObj->Delete();
 376       }
 377       else
 378       {
 379         vtkErrorMacro(<< "Unexpected error during computation of the difference.");
 380         return NULL;
 381       }
 382     }
 383     return anOutput;
 384   }
 385   else
 386   {
 387     vtkErrorMacro("We cannot yet compute difference of this type of dataset.");
 388     return NULL;
 389   }
 390 }
 391
 392 //--------------------------------------------------------------------------------------------------
 393 vtkDataSet* vtkDifferenceTimestepsFilter::DifferenceDataSet(vtkDataSet* theInput1,
 394                                                             vtkDataSet* theInput2)
 395 {
 396   vtkDataSet *anInput[2];
 397   anInput[0] = theInput1;
 398   anInput[1] = theInput2;
 399
 400   // Copy input structure into output
 401   vtkDataSet *anOutput = anInput[0]->NewInstance();
 402   anOutput->CopyStructure(anInput[0]);
 403
 404   std::vector<vtkDataArray*> anArrays;
 405   vtkDataArray *anOutputArray;
 406
 407   // Compute the difference of the the specified point or cell data array
 408   vtkDataArray* aDataArray0 = this->GetInputArrayToProcess(0, anInput[0]);
 409   vtkDataArray* aDataArray1 = this->GetInputArrayToProcess(0, anInput[1]);
 410   if (aDataArray0 == NULL || aDataArray1 == NULL)
 411   {
 412     vtkErrorMacro(<< "Input array to process is empty.");
 413     return NULL;
 414   }
 415   anArrays.push_back(aDataArray0);
 416   anArrays.push_back(aDataArray1);
 417
 418   if (anArrays.size() > 1)
 419   {
 420     if (!this->VerifyArrays(&anArrays[0], 2))
 421     {
 422       vtkErrorMacro(<< "Verification of data arrays has failed.");
 423       return NULL;
 424     }
 425
 426     anOutputArray = this->DifferenceDataArray(&anArrays[0], anArrays[0]->GetNumberOfTuples());
 427     // Determine a field association
 428     int aTypeFieldAssociation = this->GetInputFieldAssociation();
 429     if (aTypeFieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
 430     {
 431       // For point data
 432       anOutput->GetPointData()->AddArray(anOutputArray);
 433     }
 434     else if (aTypeFieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS)
 435     {
 436       // For cell data
 437       anOutput->GetCellData()->AddArray(anOutputArray);
 438     }
 439     else
 440     {
 441       vtkErrorMacro(<< "Solution is not implemeted yet.");
 442       return NULL;
 443     }
 444     anOutputArray->Delete();
 445     anArrays.clear();
 446   }
 447
 448   return anOutput;
 449 }
 450
 451 //--------------------------------------------------------------------------------------------------
 452 vtkDataArray* vtkDifferenceTimestepsFilter::DifferenceDataArray(vtkDataArray** theArrays,
 453                                                                 vtkIdType      theNumTuple)
 454 {
 455   // Create the output array based on the number of tuple and components
 456   // with a new name containing the specified prefix
 457   int aNumComp = theArrays[0]->GetNumberOfComponents();
 458   vtkDataArray *anOutput =
 459     DataTempDiffArray(theArrays[0], aNumComp, theNumTuple, this->ArrayNamePrefix);
 460
 461   // Now do the computation of the difference
 462   switch (theArrays[0]->GetDataType())
 463   {
 464     vtkTemplateMacro(
 465       vtkTemporalDataDifference(this, anOutput, theArrays, aNumComp, static_cast<VTK_TT *>(0)));
 466   default:
 467     vtkErrorMacro(<< "Execute: unknown scalar type.");
 468     return NULL;
 469   }
 470
 471   return anOutput;
 472 }
 473
 474 //--------------------------------------------------------------------------------------------------
 475 int vtkDifferenceTimestepsFilter::GetInputFieldAssociation()
 476 {
 477   vtkInformationVector *anInputArrayVec = this->GetInformation()->Get(INPUT_ARRAYS_TO_PROCESS());
 478   vtkInformation *anInputArrayInfo = anInputArrayVec->GetInformationObject(0);
 479   return anInputArrayInfo->Get(vtkDataObject::FIELD_ASSOCIATION());
 480 }
 481
 482 //--------------------------------------------------------------------------------------------------
 483 bool vtkDifferenceTimestepsFilter::VerifyArrays(vtkDataArray **theArrays, int theNumArrays)
 484 {
 485   // Get all required data to compare with other
 486   const char* anArrayName  = theArrays[0]->GetName();
 487   vtkIdType aNumTuples     = theArrays[0]->GetNumberOfTuples();
 488   vtkIdType aNumComponents = theArrays[0]->GetNumberOfComponents();
 489
 490   for (int i = 1; i < theNumArrays; ++i)
 491   {
 492     if (strcmp(theArrays[i]->GetName(), anArrayName) != 0)
 493     {
 494       vtkWarningMacro(<< "Computation of difference aborted for dataset because "
 495                       << "the array name in each time step are different.")
 496       return false;
 497     }
 498
 499     if (theArrays[i]->GetNumberOfTuples() != aNumTuples)
 500     {
 501       vtkWarningMacro(<< "Computation of difference aborted for dataset because "
 502                       << "the number of tuples in each time step are different.")
 503       return false;
 504     }
 505
 506     if (theArrays[i]->GetNumberOfComponents() != aNumComponents)
 507     {
 508       vtkWarningMacro(<< "Computation of difference aborted for dataset because "
 509                       << "the number of components in each time step are different.")
 510       return false;
 511     }
 512   }
 513
 514   return true;
 515 }