Salome HOME
Merge branch 'OCCT780'
[modules/paravis.git] / src / Plugins / MEDReader / plugin / MEDReaderIO / vtkMEDReader.cxx
1 // Copyright (C) 2010-2024  CEA, EDF
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 : Anthony Geay
20
21 #include "vtkMEDReader.h"
22 #include "vtkGenerateVectors.h"
23 #include "MEDUtilities.hxx"
24
25 #include "vtkCellArray.h"
26 #include "vtkCellData.h"
27 #include "vtkCellType.h"
28 #include "vtkDataSetAttributes.h"
29 #include "vtkDataArraySelection.h"
30 #include "vtkDoubleArray.h"
31 #include "vtkExecutive.h"
32 #include "vtkInformation.h"
33 #include "vtkInformationDataObjectMetaDataKey.h"
34 #include "vtkInformationDoubleVectorKey.h"
35 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
36 #include "vtkInformationStringKey.h"
37 #include "vtkInformationVector.h"
38 #include "vtkMultiBlockDataSet.h"
39 #include "vtkMultiTimeStepAlgorithm.h"
40 #include "vtkMutableDirectedGraph.h"
41 #include "vtkObjectFactory.h"
42 #include "vtkPointData.h"
43 #include "vtkQuadratureSchemeDefinition.h"
44 #include "vtkSmartPointer.h"
45 #include "vtkStreamingDemandDrivenPipeline.h"
46 #include "vtkStringArray.h"
47 #include "vtkUnsignedCharArray.h"
48 #include "vtkUnstructuredGrid.h"
49 #include "vtkVariantArray.h"
50
51 #ifdef MEDREADER_USE_MPI
52 #include "vtkMultiProcessController.h"
53 #include "vtkGhostCellsGenerator.h"
54 #endif
55
56 #include "MEDFileFieldRepresentationTree.hxx"
57
58 #include <map>
59 #include <string>
60 #include <vector>
61 #include <sstream>
62 #include <algorithm>
63
64 class vtkMEDReader::vtkMEDReaderInternal
65 {
66
67 public:
68   vtkMEDReaderInternal(vtkMEDReader *master):TK(0),SIL(0),LastLev0(-1)
69   {
70   }
71
72   ~vtkMEDReaderInternal()
73   {
74     if(this->SIL)
75       this->SIL->Delete();
76   }
77 public:
78   MEDFileFieldRepresentationTree Tree;
79   TimeKeeper TK;
80
81   std::string DftMeshName;
82   // Store the vtkMutableDirectedGraph that represents links between family, groups and cell types
83   vtkMutableDirectedGraph* SIL;
84   // store the lev0 id in Tree corresponding to the TIME_STEPS in the pipeline.
85   int LastLev0;
86 };
87
88 vtkStandardNewMacro(vtkMEDReader)
89
90 // vtkInformationKeyMacro(vtkMEDReader, META_DATA, DataObjectMetaData); // Here we need to customize vtkMEDReader::META_DATA method
91 // start of overload of vtkInformationKeyMacro
92 static vtkInformationDataObjectMetaDataKey *vtkMEDReader_META_DATA=new vtkInformationDataObjectMetaDataKey("META_DATA","vtkFileSeriesGroupReader");
93
94 vtkInformationDataObjectMetaDataKey *vtkMEDReader::META_DATA()
95 {
96   static const char ZE_KEY[]="vtkFileSeriesGroupReader::META_DATA";
97   vtkInformationDataObjectMetaDataKey *ret(vtkMEDReader_META_DATA);
98   MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance());
99   if(!gd->hasKey(ZE_KEY))
100     {// here META_DATA is put on global var to be exchanged with other filters without dependancy of MEDReader. Please do not change ZE_KEY !
101       std::ostringstream oss; oss << ret;
102       gd->setKeyValue(ZE_KEY,oss.str());
103     }
104   return ret;
105 }
106
107 static vtkInformationGaussDoubleVectorKey *vtkMEDReader_GAUSS_DATA=new vtkInformationGaussDoubleVectorKey("GAUSS_DATA","vtkFileSeriesGroupReader");
108
109 vtkInformationGaussDoubleVectorKey *vtkMEDReader::GAUSS_DATA()
110 {
111   static const char ZE_KEY[]="vtkFileSeriesGroupReader::GAUSS_DATA";
112   vtkInformationGaussDoubleVectorKey *ret(vtkMEDReader_GAUSS_DATA);
113   MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance());
114   if(!gd->hasKey(ZE_KEY))
115     {// here META_DATA is put on global var to be exchanged with other filters without dependancy of MEDReader. Please do not change ZE_KEY !
116       vtkInformationDoubleVectorKey *ret2(ret);
117       std::ostringstream oss; oss << ret2;
118       gd->setKeyValue(ZE_KEY,oss.str());
119     }
120   return ret;
121 }
122 // end of overload of vtkInformationKeyMacro
123
124 vtkMEDReader::vtkMEDReader():Internal(new vtkMEDReaderInternal(this))
125 {
126   this->SetNumberOfInputPorts(0);
127   this->SetNumberOfOutputPorts(1);
128 }
129
130 vtkMEDReader::~vtkMEDReader()
131 {
132   delete this->Internal;
133   this->Internal = 0;
134 }
135
136 void vtkMEDReader::Reload()
137 {
138   this->ReloadInternals();
139   this->IsStdOrMode = false;
140   this->GenerateVect = false;
141   this->GCGCP = true;
142   this->RemoveDebugArrays = false;
143   this->FieldSelection->RemoveAllArrays();
144   this->TimeFlagSelection->RemoveAllArrays();
145   this->Modified();
146 }
147 void vtkMEDReader::ReloadInternals()
148 {
149   delete this->Internal;
150   this->Internal=new vtkMEDReaderInternal(this);
151   this->Modified();
152 }
153
154 void vtkMEDReader::GenerateVectors(int val)
155 {
156   if ( !this->Internal )
157     return;
158
159   bool val2((bool)val);
160   if(val2!=this->GenerateVect)
161     {
162       this->GenerateVect=val2;
163       this->Modified();
164     }
165 }
166
167 void vtkMEDReader::ChangeMode(int newMode)
168 {
169   if ( !this->Internal )
170     return;
171
172   this->IsStdOrMode=newMode!=0;
173   this->Modified();
174 }
175
176 void vtkMEDReader::GhostCellGeneratorCallForPara(int gcgcp)
177 {
178   if ( !this->Internal )
179     return;
180
181   bool newVal(gcgcp!=0);
182   if(newVal!=this->GCGCP)
183     {
184       this->GCGCP=newVal;
185       this->Modified();
186     }
187 }
188
189 void vtkMEDReader::GetRidOffDebugArrays(int rmda)
190 {
191   if ( !this->Internal )
192     return;
193
194   bool newVal(rmda!=0);
195   if(newVal!=this->RemoveDebugArrays)
196     {
197       this->RemoveDebugArrays=newVal;
198       this->Modified();
199     }
200 }
201
202 const char *vtkMEDReader::GetSeparator()
203 {
204   return MEDFileFieldRepresentationLeavesArrays::ZE_SEP;
205 }
206
207 void vtkMEDReader::SetFileName(const char *fname)
208 {
209   if(!this->Internal)
210     return;
211   try
212     {
213       this->FileName=fname;
214       this->Modified();
215     }
216   catch(INTERP_KERNEL::Exception& e)
217     {
218       delete this->Internal;
219       this->Internal=0;
220       std::ostringstream oss;
221       oss << "Exception has been thrown in vtkMEDReader::SetFileName : " << e.what() << std::endl;
222       if(this->HasObserver("ErrorEvent") )
223         this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
224       else
225         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
226       vtkObject::BreakOnError();
227     }
228 }
229
230 char *vtkMEDReader::GetFileName()
231 {
232   if (!this->Internal)
233     return 0;
234   return const_cast<char *>(this->FileName.c_str());
235 }
236
237 int vtkMEDReader::RequestInformation(vtkInformation *request, vtkInformationVector ** /*inputVector*/, vtkInformationVector *outputVector)
238 {
239 //  std::cout << "########################################## vtkMEDReader::RequestInformation ##########################################" << std::endl;
240   if(!this->Internal)
241     return 0;
242   try
243     {
244       // Process file meta data
245       if(this->Internal->Tree.getNumberOfLeavesArrays()==0)
246         {
247           int iPart(-1),nbOfParts(-1);
248 #ifdef MEDREADER_USE_MPI
249           if (this->DistributeWithMPI)
250           {
251             vtkMultiProcessController *vmpc(vtkMultiProcessController::GetGlobalController());
252             if(vmpc)
253             {
254               iPart=vmpc->GetLocalProcessId();
255               nbOfParts=vmpc->GetNumberOfProcesses();
256             }
257           }
258 #endif
259           this->Internal->Tree.loadMainStructureOfFile(this->FileName.c_str(),iPart,nbOfParts);
260
261           // Leaves
262           this->Internal->Tree.activateTheFirst();//This line manually initialize the status of server (this) with the remote client.
263           for (int idLeaveArray = 0; idLeaveArray < this->Internal->Tree.getNumberOfLeavesArrays(); idLeaveArray++)
264           {
265             std::string name = this->Internal->Tree.getNameOf(idLeaveArray);
266             bool status = this->Internal->Tree.getStatusOf(idLeaveArray);
267             this->FieldSelection->AddArray(name.c_str(), status);
268           }
269         }
270
271       // Time flags
272       this->Internal->TK.setMaxNumberOfTimeSteps(this->Internal->Tree.getMaxNumberOfTimeSteps());
273       auto timeFlagsArray = this->Internal->TK.getTimesFlagArray();
274       for (int idTimeFlag = 0; idTimeFlag < timeFlagsArray.size() ; idTimeFlag++)
275       {
276         std::string name = timeFlagsArray[idTimeFlag].second;
277         bool status = timeFlagsArray[idTimeFlag].first;
278         this->TimeFlagSelection->AddArray(name.c_str(), status);
279       }
280
281       // Make sure internal model are synchronized
282       /// So the SIL is up to date
283       int nArrays = this->FieldSelection->GetNumberOfArrays();
284       for(int i = nArrays - 1; i >= 0; i--)
285       {
286         try
287         {
288         this->Internal->Tree.changeStatusOfAndUpdateToHaveCoherentVTKDataSet(
289           this->Internal->Tree.getIdHavingZeName(this->FieldSelection->GetArrayName(i)),
290           this->FieldSelection->GetArraySetting(i));
291         }
292         catch(INTERP_KERNEL::Exception& e)
293         {
294           // Remove the incorrect array
295           this->FieldSelection->RemoveArrayByIndex(i);
296         }
297       }
298
299 //      request->Print(cout);
300       vtkInformation *outInfo(outputVector->GetInformationObject(0));
301       outInfo->Set(vtkDataObject::DATA_TYPE_NAME(),"vtkMultiBlockDataSet");
302       this->UpdateSIL(request, outInfo);
303
304       // Set the meta data graph as a meta data key in the information
305       // That's all that is needed to transfer it along the pipeline
306       outInfo->Set(vtkMEDReader::META_DATA(),this->Internal->SIL);
307
308       bool dummy(false);
309       this->PublishTimeStepsIfNeeded(outInfo,dummy);
310     }
311   catch(INTERP_KERNEL::Exception& e)
312     {
313       std::ostringstream oss;
314       oss << "Exception has been thrown in vtkMEDReader::RequestInformation : " << e.what() << std::endl;
315       if(this->HasObserver("ErrorEvent") )
316         this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
317       else
318         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
319       vtkObject::BreakOnError();
320       return 0;
321     }
322   return 1;
323 }
324
325 int vtkMEDReader::RequestData(vtkInformation *request, vtkInformationVector ** /*inputVector*/, vtkInformationVector *outputVector)
326 {
327 //  std::cout << "########################################## vtkMEDReader::RequestData ##########################################" << std::endl;
328   if(!this->Internal)
329     return 0;
330   try
331   {
332       for(int i = 0; i < this->FieldSelection->GetNumberOfArrays(); i++)
333       {
334         this->Internal->Tree.changeStatusOfAndUpdateToHaveCoherentVTKDataSet(
335           this->Internal->Tree.getIdHavingZeName(this->FieldSelection->GetArrayName(i)),
336           this->FieldSelection->GetArraySetting(i));
337       }
338
339       auto& timeFlagsArray = this->Internal->TK.getTimesFlagArray();
340       if (timeFlagsArray.size() != this->TimeFlagSelection->GetNumberOfArrays())
341       {
342         throw INTERP_KERNEL::Exception("Unexpected size of TimeFlagSelection");
343       }
344       for(int i = 0; i < this->TimeFlagSelection->GetNumberOfArrays(); i++)
345       {
346         timeFlagsArray[i] = std::make_pair(this->TimeFlagSelection->GetArraySetting(i),
347           this->TimeFlagSelection->GetArrayName(i));
348       }
349
350 //      request->Print(cout);
351       vtkInformation *outInfo(outputVector->GetInformationObject(0));
352       vtkMultiBlockDataSet *output(vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
353       //bool isUpdated(false); // todo: unused
354       double reqTS(0.);
355       if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
356         reqTS=outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
357       ExportedTinyInfo ti;
358 #ifndef MEDREADER_USE_MPI
359       this->FillMultiBlockDataSetInstance(output,reqTS,&ti);
360 #else
361       if (this->DistributeWithMPI && this->GCGCP)
362       {
363         vtkSmartPointer<vtkGhostCellsGenerator> gcg(vtkSmartPointer<vtkGhostCellsGenerator>::New());
364         {
365           vtkDataSet *ret(RetrieveDataSetAtTime(reqTS,&ti));
366           gcg->SetInputData(ret);
367           ret->Delete();
368         }
369               // To be checked
370         // gcg->SetUseGlobalPointIds(true);
371         gcg->SetBuildIfRequired(false);
372         gcg->Update();
373         output->SetBlock(0,gcg->GetOutput());
374       }
375       else
376               this->FillMultiBlockDataSetInstance(output,reqTS,&ti);
377 #endif
378       if(!ti.empty())
379         {
380           const std::vector<double>& data(ti.getData());
381           outInfo->Set(vtkMEDReader::GAUSS_DATA(),&data[0],(int)data.size());
382           request->Append(vtkExecutive::KEYS_TO_COPY(),vtkMEDReader::GAUSS_DATA());// Thank you to SciberQuest and DIPOLE_CENTER ! Don't understand why ! In RequestInformation it does not work !
383         }
384       output->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(),reqTS);
385       // Is it really needed ? TODO
386       this->UpdateSIL(request, outInfo);
387     }
388   catch(INTERP_KERNEL::Exception& e)
389     {
390       std::cerr << "Exception has been thrown in vtkMEDReader::RequestData : " << e.what() << std::endl;
391       return 0;
392     }
393   return 1;
394 }
395
396 //------------------------------------------------------------------------------
397 int vtkMEDReader::GetNumberOfFieldsTreeArrays()
398 {
399   return this->FieldSelection->GetNumberOfArrays();
400 }
401
402 //------------------------------------------------------------------------------
403 const char* vtkMEDReader::GetFieldsTreeArrayName(int index)
404 {
405   return this->FieldSelection->GetArrayName(index);
406 }
407
408 //------------------------------------------------------------------------------
409 int vtkMEDReader::GetFieldsTreeArrayStatus(const char* name)
410 {// EDF30038 : This method can alterate this->FieldSelection !
411   return this->FieldSelection->ArrayIsEnabled(name);
412 }
413
414 //------------------------------------------------------------------------------
415 void vtkMEDReader::SetFieldsStatus(const char* name, int status)
416 {
417   // EDF30038 : GetFieldsTreeArrayStatus does not inform if the name entry already exists. So start to deal with this
418   if( ! this->FieldSelection->ArrayExists( name ) )
419   {
420     this->FieldSelection->SetArraySetting( name, status );
421     this->Modified();
422     return ;
423   }
424   if (this->GetFieldsTreeArrayStatus(name) != status)
425   {
426     if (status)
427     {
428       this->FieldSelection->EnableArray(name);
429     }
430     else
431     {
432       this->FieldSelection->DisableArray(name);
433     }
434     this->Modified();
435   }
436 }
437
438 //------------------------------------------------------------------------------
439 int vtkMEDReader::GetNumberOfTimesFlagsArrays()
440 {
441   return this->TimeFlagSelection->GetNumberOfArrays();
442 }
443
444 //------------------------------------------------------------------------------
445 const char* vtkMEDReader::GetTimesFlagsArrayName(int index)
446 {
447   return this->TimeFlagSelection->GetArrayName(index);
448 }
449
450 //------------------------------------------------------------------------------
451 int vtkMEDReader::GetTimesFlagsArrayStatus(const char* name)
452 {
453   return this->TimeFlagSelection->ArrayIsEnabled(name);
454 }
455
456 //------------------------------------------------------------------------------
457 void vtkMEDReader::SetTimesFlagsStatus(const char* name, int status)
458 {
459   if (this->GetTimesFlagsArrayStatus(name) != status)
460   {
461     if (status)
462     {
463       this->TimeFlagSelection->EnableArray(name);
464     }
465     else
466     {
467       this->TimeFlagSelection->DisableArray(name);
468     }
469     this->Modified();
470   }
471 }
472
473 void vtkMEDReader::UpdateSIL(vtkInformation* /*request*/, vtkInformation * /*info*/)
474 {
475   if(!this->Internal)
476       return;
477   std::string meshName(this->Internal->Tree.getActiveMeshName());
478   if(!this->Internal->SIL || meshName!=this->Internal->DftMeshName)
479     {
480       vtkMutableDirectedGraph *sil(vtkMutableDirectedGraph::New());
481       this->BuildSIL(sil);
482       if(this->Internal->SIL)
483         this->Internal->SIL->Delete();
484       this->Internal->SIL=sil;
485       this->Internal->DftMeshName=meshName;
486     }
487 }
488
489 /*!
490  * The returned string is the name of the mesh activated which groups and families are in \a sil.
491  */
492 std::string vtkMEDReader::BuildSIL(vtkMutableDirectedGraph* sil)
493 {
494   if (!this->Internal)
495     return std::string();
496   sil->Initialize();
497   vtkSmartPointer<vtkVariantArray> childEdge(vtkSmartPointer<vtkVariantArray>::New());
498   childEdge->InsertNextValue(0);
499   vtkSmartPointer<vtkVariantArray> crossEdge(vtkSmartPointer<vtkVariantArray>::New());
500   crossEdge->InsertNextValue(1);
501   // CrossEdge is an edge linking hierarchies.
502   vtkUnsignedCharArray* crossEdgesArray=vtkUnsignedCharArray::New();
503   crossEdgesArray->SetName("CrossEdges");
504   sil->GetEdgeData()->AddArray(crossEdgesArray);
505   crossEdgesArray->Delete();
506   std::vector<std::string> names;
507   // Now build the hierarchy.
508   vtkIdType rootId=sil->AddVertex();
509   names.push_back("SIL");
510   // Add global fields root
511   vtkIdType fieldsRoot(sil->AddChild(rootId,childEdge));
512   names.push_back("FieldsStatusTree");
513   this->Internal->Tree.feedSIL(sil,fieldsRoot,childEdge,names);
514   vtkIdType meshesFamsGrpsRoot(sil->AddChild(rootId,childEdge));
515   names.push_back("MeshesFamsGrps");
516   std::string dftMeshName(this->Internal->Tree.feedSILForFamsAndGrps(sil,meshesFamsGrpsRoot,childEdge,names));
517   // This array is used to assign names to nodes.
518   vtkStringArray *namesArray(vtkStringArray::New());
519   namesArray->SetName("Names");
520   namesArray->SetNumberOfTuples(sil->GetNumberOfVertices());
521   sil->GetVertexData()->AddArray(namesArray);
522   namesArray->Delete();
523   std::vector<std::string>::const_iterator iter;
524   vtkIdType cc;
525   for(cc=0, iter=names.begin(); iter!=names.end(); ++iter, ++cc)
526     namesArray->SetValue(cc,(*iter).c_str());
527   return dftMeshName;
528 }
529
530 double vtkMEDReader::PublishTimeStepsIfNeeded(vtkInformation *outInfo, bool& isUpdated)
531 {
532   if(!this->Internal)
533     return 0.0;
534
535   int lev0(-1);
536   std::vector<double> tsteps;
537   if(!this->IsStdOrMode)
538     tsteps=this->Internal->Tree.getTimeSteps(lev0,this->Internal->TK);
539   else
540     { tsteps.resize(1); tsteps[0]=0.; }
541   isUpdated=false;
542   if(lev0!=this->Internal->LastLev0)
543     {
544       isUpdated=true;
545       double timeRange[2];
546       timeRange[0]=tsteps.front();
547       timeRange[1]=tsteps.back();
548       outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),&tsteps[0],(int)tsteps.size());
549       outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),timeRange,2);
550       this->Internal->LastLev0=lev0;
551     }
552   return tsteps.front();
553 }
554
555 void vtkMEDReader::FillMultiBlockDataSetInstance(vtkMultiBlockDataSet *output, double reqTS, ExportedTinyInfo *internalInfo)
556 {
557   if( !this->Internal )
558     return;
559   vtkDataSet *ret(RetrieveDataSetAtTime(reqTS,internalInfo));
560   output->SetBlock(0,ret);
561   ret->Delete();
562 }
563
564 vtkDataSet *vtkMEDReader::RetrieveDataSetAtTime(double reqTS, ExportedTinyInfo *internalInfo)
565 {
566   if( !this->Internal )
567     return 0;
568   std::string meshName;
569   vtkDataSet *ret(this->Internal->Tree.buildVTKInstance(this->IsStdOrMode,reqTS,meshName,this->Internal->TK,!this->RemoveDebugArrays,internalInfo));
570   if(this->GenerateVect)
571     {
572       vtkGenerateVectors::Operate(ret->GetPointData());
573       vtkGenerateVectors::Operate(ret->GetCellData());
574       vtkGenerateVectors::Operate(ret->GetFieldData());
575       // The operations above have potentially created new arrays -> This breaks the optimization of StaticMesh that expects the same field arrays over time.
576       // To enforce the cache recomputation declare modification of mesh.
577       //vtkGenerateVectors::ChangeMeshTimeToUpdateCache(ret);
578     }
579   return ret;
580 }
581
582 void vtkMEDReader::PrintSelf(ostream& os, vtkIndent indent)
583 {
584   this->Superclass::PrintSelf(os, indent);
585 }