1 // Copyright (C) 2014-2023 CEA, EDF
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Maxim Glibin
21 #include "vtkDifferenceTimestepsFilter.h"
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>
37 // Temporal difference of data array
38 vtkDataArray* DataTempDiffArray(
39 vtkDataArray* theDataArray, vtkIdType theNumComp, vtkIdType theNumTuple, const char* thePrefix)
41 // Create the new array
42 vtkAbstractArray* anAbstractArray = theDataArray->CreateArray(theDataArray->GetDataType());
43 vtkDataArray* anOutput = vtkDataArray::SafeDownCast(anAbstractArray);
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());
54 // Templated difference function
56 void vtkTemporalDataDifference(vtkDifferenceTimestepsFilter* /*theDTF*/, vtkDataArray* theOutput,
57 vtkDataArray** theArrays, vtkIdType theNumComp, T*)
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));
63 vtkIdType N = theArrays[0]->GetNumberOfTuples();
64 for (vtkIdType t = 0; t < N; ++t)
66 T* x0 = &anInputData0[t * theNumComp];
67 T* x1 = &anInputData1[t * theNumComp];
68 for (int c = 0; c < theNumComp; ++c)
70 // Compute the difference
71 *anOutputData++ = static_cast<T>(x1[c] - x0[c]);
74 // Copy component name
75 for (int c = 0; c < theNumComp; ++c)
77 theOutput->SetComponentName(c, theArrays[0]->GetComponentName(c));
79 theOutput->SetNumberOfTuples(N);
82 vtkStandardNewMacro(vtkDifferenceTimestepsFilter)
84 //--------------------------------------------------------------------------------------------------
85 vtkDifferenceTimestepsFilter::vtkDifferenceTimestepsFilter()
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;
95 this->SetNumberOfInputPorts(1);
96 this->SetNumberOfOutputPorts(1);
98 // Set the input data array that the algorithm will process
99 this->SetInputArrayToProcess(
100 0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
103 //--------------------------------------------------------------------------------------------------
104 vtkDifferenceTimestepsFilter::~vtkDifferenceTimestepsFilter()
106 this->TimeStepValues.clear();
107 this->SetArrayNamePrefix(nullptr);
110 //--------------------------------------------------------------------------------------------------
111 void vtkDifferenceTimestepsFilter::PrintSelf(ostream& theOS, vtkIndent theIndent)
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;
121 //--------------------------------------------------------------------------------------------------
122 int vtkDifferenceTimestepsFilter::FillInputPortInformation(int thePort, vtkInformation* theInfo)
125 theInfo->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
130 //--------------------------------------------------------------------------------------------------
131 int vtkDifferenceTimestepsFilter::FillOutputPortInformation(
132 int vtkNotUsed(thePort), vtkInformation* theInfo)
134 theInfo->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkDataObject");
138 //--------------------------------------------------------------------------------------------------
139 int vtkDifferenceTimestepsFilter::RequestDataObject(vtkInformation* vtkNotUsed(theRequest),
140 vtkInformationVector** theInputVector, vtkInformationVector* theOutputVector)
142 if (this->GetNumberOfInputPorts() == 0 || this->GetNumberOfOutputPorts() == 0)
145 vtkInformation* anInputInfo = theInputVector[0]->GetInformationObject(0);
146 if (anInputInfo == nullptr)
148 vtkErrorMacro(<< "Input information vector is missed.");
152 vtkDataObject* anInputObj = anInputInfo->Get(vtkDataObject::DATA_OBJECT());
153 if (anInputObj != nullptr)
156 for (int i = 0; i < this->GetNumberOfOutputPorts(); ++i)
158 vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(i);
159 vtkDataObject* anOutputObj = anOutputInfo->Get(vtkDataObject::DATA_OBJECT());
160 if (!anOutputObj || !anOutputObj->IsA(anInputObj->GetClassName()))
162 vtkDataObject* aNewOutput = anInputObj->NewInstance();
163 anOutputInfo->Set(vtkDataObject::DATA_OBJECT(), aNewOutput);
164 aNewOutput->Delete();
172 //--------------------------------------------------------------------------------------------------
173 int vtkDifferenceTimestepsFilter::RequestInformation(vtkInformation* vtkNotUsed(theRequest),
174 vtkInformationVector** theInputVector, vtkInformationVector* /*theOutputVector*/)
176 // Get input and output information objects
177 vtkInformation* anInInfo = theInputVector[0]->GetInformationObject(0);
178 //vtkInformation* anOutInfo = theOutputVector->GetInformationObject(0); // todo: unused
180 // Check for presence more than one time step
181 if (anInInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
183 // Find time on input
184 this->NumberTimeSteps = anInInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
185 if (this->NumberTimeSteps < 2)
187 vtkErrorMacro(<< "Not enough numbers of time steps: " << this->NumberTimeSteps);
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)
195 vtkErrorMacro(<< "Array of time steps is empty.");
201 vtkErrorMacro(<< "No time steps in input data.");
205 // Update range of indices of the time steps
206 this->RangeIndicesTimeSteps[0] = 0;
207 this->RangeIndicesTimeSteps[1] = this->NumberTimeSteps - 1;
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 ...
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());
222 // Unset the time range
223 if(anOutInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_RANGE()))
224 anOutInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
229 //--------------------------------------------------------------------------------------------------
230 int vtkDifferenceTimestepsFilter::RequestUpdateExtent(vtkInformation* /*theRequest*/,
231 vtkInformationVector** theInputVector, vtkInformationVector* theOutputVector)
233 // Get the information objects
234 vtkInformation* anInputInfo = theInputVector[0]->GetInformationObject(0);
235 vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(0);
237 // Indices must not go beyond the range of indices of the time steps
238 if (this->FirstTimeStepIndex >= this->NumberTimeSteps || this->FirstTimeStepIndex < 0)
240 vtkErrorMacro(<< "Specified index of the first time step [" << this->FirstTimeStepIndex
241 << "] is outside the range of indices.");
244 if (this->SecondTimeStepIndex >= this->NumberTimeSteps || this->SecondTimeStepIndex < 0)
246 vtkErrorMacro(<< "Specified index of the second time step [" << this->SecondTimeStepIndex
247 << "] is outside the range of indices.");
251 // Warn if the selected time steps are equal
252 if (this->FirstTimeStepIndex == this->SecondTimeStepIndex)
254 vtkWarningMacro(<< "First and second indices [" << this->FirstTimeStepIndex << " = "
255 << this->SecondTimeStepIndex << "] are the same.");
258 // Find the required input time steps and request them
259 if (anOutputInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
261 // Get the available input times
262 double* anInputTimes = anInputInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
263 if (anInputTimes != nullptr)
265 // Compute the requested times
266 double anInputUpdateTimes[2];
267 int aNumberInputUpdateTimes(0);
269 // For each the requested time mark the required input times
270 anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->FirstTimeStepIndex];
271 anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->SecondTimeStepIndex];
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);
282 //--------------------------------------------------------------------------------------------------
283 /*int vtkDifferenceTimestepsFilter::RequestData(vtkInformation* vtkNotUsed(theRequest),
284 vtkInformationVector** theInputVector, vtkInformationVector* theOutputVector)
286 // Get the information objects
287 vtkInformation* anInputInfo = theInputVector[0]->GetInformationObject(0);
288 vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(0);
290 vtkDataObject* anOutputDataObj = nullptr;
292 vtkMultiBlockDataSet* anInputData =
293 vtkMultiBlockDataSet::SafeDownCast(anInputInfo->Get(vtkDataObject::DATA_OBJECT()));
295 int aNumberTimeSteps = anInputData->GetNumberOfBlocks();
296 if (aNumberTimeSteps == 2)
299 vtkDataObject* aData0 = anInputData->GetBlock(0);
300 vtkDataObject* aData1 = anInputData->GetBlock(1);
301 if (aData0 == nullptr && aData1 == nullptr)
303 vtkErrorMacro(<< "Null data set.");
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();
315 vtkErrorMacro(<< "The amount of time blocks is not correct: " << aNumberTimeSteps);
322 //--------------------------------------------------------------------------------------------------
323 vtkDataObject* vtkDifferenceTimestepsFilter::DifferenceDataObject(
324 vtkDataObject* theInput1, vtkDataObject* theInput2)
326 // Determine the input object type
327 if (theInput1->IsA("vtkDataSet"))
329 vtkDataSet* anInDataSet1 = vtkDataSet::SafeDownCast(theInput1);
330 vtkDataSet* anInDataSet2 = vtkDataSet::SafeDownCast(theInput2);
331 return this->DifferenceDataSet(anInDataSet1, anInDataSet2);
333 else if (theInput1->IsA("vtkCompositeDataSet"))
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);
340 vtkCompositeDataSet* anOutput = aMGDataSet[0]->NewInstance();
341 anOutput->CopyStructure(aMGDataSet[0]);
343 vtkSmartPointer<vtkCompositeDataIterator> anIter;
344 anIter.TakeReference(aMGDataSet[0]->NewIterator());
345 for (anIter->InitTraversal(); !anIter->IsDoneWithTraversal(); anIter->GoToNextItem())
347 vtkDataObject* aDataObj1 = anIter->GetCurrentDataObject();
348 vtkDataObject* aDataObj2 = aMGDataSet[1]->GetDataSet(anIter);
349 if (aDataObj1 == nullptr || aDataObj2 == nullptr)
351 vtkWarningMacro("The composite datasets were not identical in structure.");
355 vtkDataObject* aResultDObj = this->DifferenceDataObject(aDataObj1, aDataObj2);
356 if (aResultDObj != nullptr)
358 anOutput->SetDataSet(anIter, aResultDObj);
359 aResultDObj->Delete();
363 vtkErrorMacro(<< "Unexpected error during computation of the difference.");
371 vtkErrorMacro("We cannot yet compute difference of this type of dataset.");
376 //--------------------------------------------------------------------------------------------------
377 vtkDataSet* vtkDifferenceTimestepsFilter::DifferenceDataSet(
378 vtkDataSet* theInput1, vtkDataSet* theInput2)
380 vtkDataSet* anInput[2];
381 anInput[0] = theInput1;
382 anInput[1] = theInput2;
384 // Copy input structure into output
385 vtkDataSet* anOutput = anInput[0]->NewInstance();
386 anOutput->CopyStructure(anInput[0]);
388 std::vector<vtkDataArray*> anArrays;
389 vtkDataArray* anOutputArray;
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)
396 vtkErrorMacro(<< "Input array to process is empty.");
399 anArrays.push_back(aDataArray0);
400 anArrays.push_back(aDataArray1);
402 if (anArrays.size() > 1)
404 if (!this->VerifyArrays(&anArrays[0], 2))
406 vtkErrorMacro(<< "Verification of data arrays has failed.");
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)
416 anOutput->GetPointData()->AddArray(anOutputArray);
418 else if (aTypeFieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS)
421 anOutput->GetCellData()->AddArray(anOutputArray);
425 vtkErrorMacro(<< "Solution is not implemeted yet.");
428 anOutputArray->Delete();
435 //--------------------------------------------------------------------------------------------------
436 vtkDataArray* vtkDifferenceTimestepsFilter::DifferenceDataArray(
437 vtkDataArray** theArrays, vtkIdType theNumTuple)
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);
445 // Now do the computation of the difference
446 switch (theArrays[0]->GetDataType())
449 vtkTemporalDataDifference(this, anOutput, theArrays, aNumComp, static_cast<VTK_TT*>(0)));
451 vtkErrorMacro(<< "Execute: unknown scalar type.");
458 //--------------------------------------------------------------------------------------------------
459 int vtkDifferenceTimestepsFilter::GetInputFieldAssociation()
461 vtkInformationVector* anInputArrayVec = this->GetInformation()->Get(INPUT_ARRAYS_TO_PROCESS());
462 vtkInformation* anInputArrayInfo = anInputArrayVec->GetInformationObject(0);
463 return anInputArrayInfo->Get(vtkDataObject::FIELD_ASSOCIATION());
466 //--------------------------------------------------------------------------------------------------
467 bool vtkDifferenceTimestepsFilter::VerifyArrays(vtkDataArray** theArrays, int theNumArrays)
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();
474 for (int i = 1; i < theNumArrays; ++i)
476 if (strcmp(theArrays[i]->GetName(), anArrayName) != 0)
478 vtkWarningMacro(<< "Computation of difference aborted for dataset because "
479 << "the array name in each time step are different.");
483 if (theArrays[i]->GetNumberOfTuples() != aNumTuples)
485 vtkWarningMacro(<< "Computation of difference aborted for dataset because "
486 << "the number of tuples in each time step are different.");
490 if (theArrays[i]->GetNumberOfComponents() != aNumComponents)
492 vtkWarningMacro(<< "Computation of difference aborted for dataset because "
493 << "the number of components in each time step are different.");