src/Plugins/DifferenceTimesteps/plugin/DifferenceTimestepsModule/vtkDifferenceTimestepsFilter.cxx

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