1 // Copyright (C) 2014-2016 CEA/DEN, EDF R&D
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(vtkDataArray* theDataArray,
40 vtkIdType theNumTuple,
41 const char* thePrefix)
43 // Create the new array
44 vtkAbstractArray *anAbstractArray = theDataArray->CreateArray(theDataArray->GetDataType());
45 vtkDataArray *anOutput = vtkDataArray::SafeDownCast(anAbstractArray);
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());
56 // Templated difference function
58 void vtkTemporalDataDifference(
59 vtkDifferenceTimestepsFilter* theDTF,
60 vtkDataArray* theOutput,
61 vtkDataArray** theArrays,
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));
69 vtkIdType N = theArrays[0]->GetNumberOfTuples();
70 for (vtkIdType t = 0; t < N; ++t)
72 T* x0 = &anInputData0[t*theNumComp];
73 T* x1 = &anInputData1[t*theNumComp];
74 for (int c = 0; c < theNumComp; ++c)
76 // Compute the difference
77 *anOutputData++ = static_cast<T>(x1[c]-x0[c]);
80 theOutput->SetNumberOfTuples(N);
83 vtkStandardNewMacro(vtkDifferenceTimestepsFilter);
85 //--------------------------------------------------------------------------------------------------
86 vtkDifferenceTimestepsFilter::vtkDifferenceTimestepsFilter()
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;
96 this->SetNumberOfInputPorts(1);
97 this->SetNumberOfOutputPorts(1);
99 // Set the input data array that the algorithm will process
100 this->SetInputArrayToProcess(
104 vtkDataObject::FIELD_ASSOCIATION_POINTS,
105 vtkDataSetAttributes::SCALARS);
108 //--------------------------------------------------------------------------------------------------
109 vtkDifferenceTimestepsFilter::~vtkDifferenceTimestepsFilter()
111 this->TimeStepValues.clear();
112 this->SetArrayNamePrefix(NULL);
115 //--------------------------------------------------------------------------------------------------
116 void vtkDifferenceTimestepsFilter::PrintSelf(ostream& theOS, vtkIndent theIndent)
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;
126 //--------------------------------------------------------------------------------------------------
127 int vtkDifferenceTimestepsFilter::FillInputPortInformation(
128 int thePort, vtkInformation* theInfo)
131 theInfo->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
136 //--------------------------------------------------------------------------------------------------
137 int vtkDifferenceTimestepsFilter::FillOutputPortInformation(
138 int vtkNotUsed(thePort), vtkInformation* theInfo)
140 theInfo->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkDataObject");
144 //--------------------------------------------------------------------------------------------------
145 int vtkDifferenceTimestepsFilter::RequestDataObject(
146 vtkInformation* vtkNotUsed(theRequest),
147 vtkInformationVector** theInputVector,
148 vtkInformationVector* theOutputVector)
150 if (this->GetNumberOfInputPorts() == 0 || this->GetNumberOfOutputPorts() == 0)
153 vtkInformation* anInputInfo = theInputVector[0]->GetInformationObject(0);
154 if (anInputInfo == NULL)
156 vtkErrorMacro(<< "Input information vector is missed.");
160 vtkDataObject *anInputObj = anInputInfo->Get(vtkDataObject::DATA_OBJECT());
161 if (anInputObj != NULL)
164 for (int i = 0; i < this->GetNumberOfOutputPorts(); ++i)
166 vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(i);
167 vtkDataObject *anOutputObj = anOutputInfo->Get(vtkDataObject::DATA_OBJECT());
168 if (!anOutputObj || !anOutputObj->IsA(anInputObj->GetClassName()))
170 vtkDataObject* aNewOutput = anInputObj->NewInstance();
171 anOutputInfo->Set(vtkDataObject::DATA_OBJECT(), aNewOutput);
172 aNewOutput->Delete();
180 //--------------------------------------------------------------------------------------------------
181 int vtkDifferenceTimestepsFilter::RequestInformation(
182 vtkInformation* vtkNotUsed(theRequest),
183 vtkInformationVector** theInputVector,
184 vtkInformationVector* theOutputVector)
186 // Get input and output information objects
187 vtkInformation *anInInfo = theInputVector[0]->GetInformationObject(0);
188 vtkInformation *anOutInfo = theOutputVector->GetInformationObject(0);
190 // Check for presence more than one time step
191 if (anInInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
193 // Find time on input
194 this->NumberTimeSteps = anInInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
195 if (this->NumberTimeSteps < 2)
197 vtkErrorMacro(<< "Not enough numbers of time steps: " << this->NumberTimeSteps);
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)
205 vtkErrorMacro(<<"Array of time steps is empty.");
211 vtkErrorMacro(<< "No time steps in input data.");
215 // Update range of indices of the time steps
216 this->RangeIndicesTimeSteps[0] = 0;
217 this->RangeIndicesTimeSteps[1] = this->NumberTimeSteps - 1;
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 ...
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());
232 // Unset the time range
233 if(anOutInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_RANGE()))
234 anOutInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
239 //--------------------------------------------------------------------------------------------------
240 int vtkDifferenceTimestepsFilter::RequestUpdateExtent(vtkInformation* theRequest,
241 vtkInformationVector** theInputVector,
242 vtkInformationVector* theOutputVector)
244 // Get the information objects
245 vtkInformation* anInputInfo = theInputVector[0]->GetInformationObject(0);
246 vtkInformation* anOutputInfo = theOutputVector->GetInformationObject(0);
248 // Indices must not go beyond the range of indices of the time steps
249 if (this->FirstTimeStepIndex >= this->NumberTimeSteps || this->FirstTimeStepIndex < 0)
251 vtkErrorMacro(<< "Specified index of the first time step ["
252 << this->FirstTimeStepIndex
253 << "] is outside the range of indices.");
256 if (this->SecondTimeStepIndex >= this->NumberTimeSteps || this->SecondTimeStepIndex < 0)
258 vtkErrorMacro(<< "Specified index of the second time step ["
259 << this->SecondTimeStepIndex
260 << "] is outside the range of indices.");
264 // Warn if the selected time steps are equal
265 if (this->FirstTimeStepIndex == this->SecondTimeStepIndex)
267 vtkWarningMacro(<< "First and second indices ["
268 << this->FirstTimeStepIndex
269 << " = " << this->SecondTimeStepIndex
270 << "] are the same.");
273 // Find the required input time steps and request them
274 if (anOutputInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
276 // Get the available input times
277 double *anInputTimes = anInputInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
278 if (anInputTimes != NULL)
280 // Compute the requested times
281 double anInputUpdateTimes[2];
282 int aNumberInputUpdateTimes(0);
284 // For each the requested time mark the required input times
285 anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->FirstTimeStepIndex];
286 anInputUpdateTimes[aNumberInputUpdateTimes++] = anInputTimes[this->SecondTimeStepIndex];
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);
297 //--------------------------------------------------------------------------------------------------
298 int vtkDifferenceTimestepsFilter::RequestData(vtkInformation* vtkNotUsed(theRequest),
299 vtkInformationVector** theInputVector,
300 vtkInformationVector* theOutputVector)
302 // Get the information objects
303 vtkInformation *anInputInfo = theInputVector[0]->GetInformationObject(0);
304 vtkInformation *anOutputInfo = theOutputVector->GetInformationObject(0);
306 vtkDataObject *anOutputDataObj = NULL;
308 vtkMultiBlockDataSet *anInputData =
309 vtkMultiBlockDataSet::SafeDownCast(anInputInfo->Get(vtkDataObject::DATA_OBJECT()));
311 int aNumberTimeSteps = anInputData->GetNumberOfBlocks();
312 if (aNumberTimeSteps == 2)
315 vtkDataObject* aData0 = anInputData->GetBlock(0);
316 vtkDataObject* aData1 = anInputData->GetBlock(1);
317 if (aData0 == NULL && aData1 == NULL)
319 vtkErrorMacro(<< "Null data set.");
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();
331 vtkErrorMacro(<< "The amount of time blocks is not correct: " << aNumberTimeSteps);
338 //--------------------------------------------------------------------------------------------------
339 vtkDataObject* vtkDifferenceTimestepsFilter::DifferenceDataObject(vtkDataObject* theInput1,
340 vtkDataObject* theInput2)
342 // Determine the input object type
343 if (theInput1->IsA("vtkDataSet"))
345 vtkDataSet *anInDataSet1 = vtkDataSet::SafeDownCast(theInput1);
346 vtkDataSet *anInDataSet2 = vtkDataSet::SafeDownCast(theInput2);
347 return this->DifferenceDataSet(anInDataSet1, anInDataSet2);
349 else if (theInput1->IsA("vtkCompositeDataSet"))
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);
356 vtkCompositeDataSet *anOutput = aMGDataSet[0]->NewInstance();
357 anOutput->CopyStructure(aMGDataSet[0]);
359 vtkSmartPointer<vtkCompositeDataIterator> anIter;
360 anIter.TakeReference(aMGDataSet[0]->NewIterator());
361 for (anIter->InitTraversal(); !anIter->IsDoneWithTraversal(); anIter->GoToNextItem())
363 vtkDataObject* aDataObj1 = anIter->GetCurrentDataObject();
364 vtkDataObject* aDataObj2 = aMGDataSet[1]->GetDataSet(anIter);
365 if (aDataObj1 == NULL || aDataObj2 == NULL)
367 vtkWarningMacro("The composite datasets were not identical in structure.");
371 vtkDataObject *aResultDObj = this->DifferenceDataObject(aDataObj1, aDataObj2);
372 if (aResultDObj != NULL)
374 anOutput->SetDataSet(anIter, aResultDObj);
375 aResultDObj->Delete();
379 vtkErrorMacro(<< "Unexpected error during computation of the difference.");
387 vtkErrorMacro("We cannot yet compute difference of this type of dataset.");
392 //--------------------------------------------------------------------------------------------------
393 vtkDataSet* vtkDifferenceTimestepsFilter::DifferenceDataSet(vtkDataSet* theInput1,
394 vtkDataSet* theInput2)
396 vtkDataSet *anInput[2];
397 anInput[0] = theInput1;
398 anInput[1] = theInput2;
400 // Copy input structure into output
401 vtkDataSet *anOutput = anInput[0]->NewInstance();
402 anOutput->CopyStructure(anInput[0]);
404 vtkstd::vector<vtkDataArray*> anArrays;
405 vtkDataArray *anOutputArray;
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)
412 vtkErrorMacro(<< "Input array to process is empty.");
415 anArrays.push_back(aDataArray0);
416 anArrays.push_back(aDataArray1);
418 if (anArrays.size() > 1)
420 if (!this->VerifyArrays(&anArrays[0], 2))
422 vtkErrorMacro(<< "Verification of data arrays has failed.");
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)
432 anOutput->GetPointData()->AddArray(anOutputArray);
434 else if (aTypeFieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS)
437 anOutput->GetCellData()->AddArray(anOutputArray);
441 vtkErrorMacro(<< "Solution is not implemeted yet.");
444 anOutputArray->Delete();
451 //--------------------------------------------------------------------------------------------------
452 vtkDataArray* vtkDifferenceTimestepsFilter::DifferenceDataArray(vtkDataArray** theArrays,
453 vtkIdType theNumTuple)
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);
461 // Now do the computation of the difference
462 switch (theArrays[0]->GetDataType())
465 vtkTemporalDataDifference(this, anOutput, theArrays, aNumComp, static_cast<VTK_TT *>(0)));
467 vtkErrorMacro(<< "Execute: unknown scalar type.");
474 //--------------------------------------------------------------------------------------------------
475 int vtkDifferenceTimestepsFilter::GetInputFieldAssociation()
477 vtkInformationVector *anInputArrayVec = this->GetInformation()->Get(INPUT_ARRAYS_TO_PROCESS());
478 vtkInformation *anInputArrayInfo = anInputArrayVec->GetInformationObject(0);
479 return anInputArrayInfo->Get(vtkDataObject::FIELD_ASSOCIATION());
482 //--------------------------------------------------------------------------------------------------
483 bool vtkDifferenceTimestepsFilter::VerifyArrays(vtkDataArray **theArrays, int theNumArrays)
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();
490 for (int i = 1; i < theNumArrays; ++i)
492 if (strcmp(theArrays[i]->GetName(), anArrayName) != 0)
494 vtkWarningMacro(<< "Computation of difference aborted for dataset because "
495 << "the array name in each time step are different.")
499 if (theArrays[i]->GetNumberOfTuples() != aNumTuples)
501 vtkWarningMacro(<< "Computation of difference aborted for dataset because "
502 << "the number of tuples in each time step are different.")
506 if (theArrays[i]->GetNumberOfComponents() != aNumComponents)
508 vtkWarningMacro(<< "Computation of difference aborted for dataset because "
509 << "the number of components in each time step are different.")