1 // Copyright (C) 2010-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 : Anthony Geay
21 #include "vtkMEDReader.h"
22 #include "vtkGenerateVectors.h"
23 #include "MEDUtilities.hxx"
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"
51 #ifdef MEDREADER_USE_MPI
52 #include "vtkMultiProcessController.h"
53 #include "vtkGhostCellsGenerator.h"
56 #include "MEDFileFieldRepresentationTree.hxx"
64 class vtkMEDReader::vtkMEDReaderInternal
68 vtkMEDReaderInternal(vtkMEDReader *master):TK(0),SIL(0),LastLev0(-1)
72 ~vtkMEDReaderInternal()
78 MEDFileFieldRepresentationTree Tree;
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.
88 vtkStandardNewMacro(vtkMEDReader)
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");
94 vtkInformationDataObjectMetaDataKey *vtkMEDReader::META_DATA()
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());
107 static vtkInformationGaussDoubleVectorKey *vtkMEDReader_GAUSS_DATA=new vtkInformationGaussDoubleVectorKey("GAUSS_DATA","vtkFileSeriesGroupReader");
109 vtkInformationGaussDoubleVectorKey *vtkMEDReader::GAUSS_DATA()
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());
122 // end of overload of vtkInformationKeyMacro
124 vtkMEDReader::vtkMEDReader():Internal(new vtkMEDReaderInternal(this))
126 this->SetNumberOfInputPorts(0);
127 this->SetNumberOfOutputPorts(1);
130 vtkMEDReader::~vtkMEDReader()
132 delete this->Internal;
136 void vtkMEDReader::Reload()
138 this->ReloadInternals();
139 this->IsStdOrMode = false;
140 this->GenerateVect = false;
142 this->RemoveDebugArrays = false;
143 this->FieldSelection->RemoveAllArrays();
144 this->TimeFlagSelection->RemoveAllArrays();
147 void vtkMEDReader::ReloadInternals()
149 delete this->Internal;
150 this->Internal=new vtkMEDReaderInternal(this);
154 void vtkMEDReader::GenerateVectors(int val)
156 if ( !this->Internal )
159 bool val2((bool)val);
160 if(val2!=this->GenerateVect)
162 this->GenerateVect=val2;
167 void vtkMEDReader::ChangeMode(int newMode)
169 if ( !this->Internal )
172 this->IsStdOrMode=newMode!=0;
176 void vtkMEDReader::GhostCellGeneratorCallForPara(int gcgcp)
178 if ( !this->Internal )
181 bool newVal(gcgcp!=0);
182 if(newVal!=this->GCGCP)
189 void vtkMEDReader::GetRidOffDebugArrays(int rmda)
191 if ( !this->Internal )
194 bool newVal(rmda!=0);
195 if(newVal!=this->RemoveDebugArrays)
197 this->RemoveDebugArrays=newVal;
202 const char *vtkMEDReader::GetSeparator()
204 return MEDFileFieldRepresentationLeavesArrays::ZE_SEP;
207 void vtkMEDReader::SetFileName(const char *fname)
213 this->FileName=fname;
216 catch(INTERP_KERNEL::Exception& e)
218 delete this->Internal;
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()));
225 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
226 vtkObject::BreakOnError();
230 char *vtkMEDReader::GetFileName()
234 return const_cast<char *>(this->FileName.c_str());
237 int vtkMEDReader::RequestInformation(vtkInformation *request, vtkInformationVector ** /*inputVector*/, vtkInformationVector *outputVector)
239 // std::cout << "########################################## vtkMEDReader::RequestInformation ##########################################" << std::endl;
244 // Process file meta data
245 if(this->Internal->Tree.getNumberOfLeavesArrays()==0)
247 int iPart(-1),nbOfParts(-1);
248 #ifdef MEDREADER_USE_MPI
249 if (this->DistributeWithMPI)
251 vtkMultiProcessController *vmpc(vtkMultiProcessController::GetGlobalController());
254 iPart=vmpc->GetLocalProcessId();
255 nbOfParts=vmpc->GetNumberOfProcesses();
259 this->Internal->Tree.loadMainStructureOfFile(this->FileName.c_str(),iPart,nbOfParts);
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++)
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);
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++)
276 std::string name = timeFlagsArray[idTimeFlag].second;
277 bool status = timeFlagsArray[idTimeFlag].first;
278 this->TimeFlagSelection->AddArray(name.c_str(), status);
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--)
288 this->Internal->Tree.changeStatusOfAndUpdateToHaveCoherentVTKDataSet(
289 this->Internal->Tree.getIdHavingZeName(this->FieldSelection->GetArrayName(i)),
290 this->FieldSelection->GetArraySetting(i));
292 catch(INTERP_KERNEL::Exception& e)
294 // Remove the incorrect array
295 this->FieldSelection->RemoveArrayByIndex(i);
299 // request->Print(cout);
300 vtkInformation *outInfo(outputVector->GetInformationObject(0));
301 outInfo->Set(vtkDataObject::DATA_TYPE_NAME(),"vtkMultiBlockDataSet");
302 this->UpdateSIL(request, outInfo);
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);
309 this->PublishTimeStepsIfNeeded(outInfo,dummy);
311 catch(INTERP_KERNEL::Exception& e)
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()));
318 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
319 vtkObject::BreakOnError();
325 int vtkMEDReader::RequestData(vtkInformation *request, vtkInformationVector ** /*inputVector*/, vtkInformationVector *outputVector)
327 // std::cout << "########################################## vtkMEDReader::RequestData ##########################################" << std::endl;
332 for(int i = 0; i < this->FieldSelection->GetNumberOfArrays(); i++)
334 this->Internal->Tree.changeStatusOfAndUpdateToHaveCoherentVTKDataSet(
335 this->Internal->Tree.getIdHavingZeName(this->FieldSelection->GetArrayName(i)),
336 this->FieldSelection->GetArraySetting(i));
339 auto& timeFlagsArray = this->Internal->TK.getTimesFlagArray();
340 if (timeFlagsArray.size() != this->TimeFlagSelection->GetNumberOfArrays())
342 throw INTERP_KERNEL::Exception("Unexpected size of TimeFlagSelection");
344 for(int i = 0; i < this->TimeFlagSelection->GetNumberOfArrays(); i++)
346 timeFlagsArray[i] = std::make_pair(this->TimeFlagSelection->GetArraySetting(i),
347 this->TimeFlagSelection->GetArrayName(i));
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
355 if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
356 reqTS=outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
358 #ifndef MEDREADER_USE_MPI
359 this->FillMultiBlockDataSetInstance(output,reqTS,&ti);
361 if (this->DistributeWithMPI && this->GCGCP)
363 vtkSmartPointer<vtkGhostCellsGenerator> gcg(vtkSmartPointer<vtkGhostCellsGenerator>::New());
365 vtkDataSet *ret(RetrieveDataSetAtTime(reqTS,&ti));
366 gcg->SetInputData(ret);
370 // gcg->SetUseGlobalPointIds(true);
371 gcg->SetBuildIfRequired(false);
373 output->SetBlock(0,gcg->GetOutput());
376 this->FillMultiBlockDataSetInstance(output,reqTS,&ti);
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 !
384 output->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(),reqTS);
385 // Is it really needed ? TODO
386 this->UpdateSIL(request, outInfo);
388 catch(INTERP_KERNEL::Exception& e)
390 std::cerr << "Exception has been thrown in vtkMEDReader::RequestData : " << e.what() << std::endl;
396 //------------------------------------------------------------------------------
397 int vtkMEDReader::GetNumberOfFieldsTreeArrays()
399 return this->FieldSelection->GetNumberOfArrays();
402 //------------------------------------------------------------------------------
403 const char* vtkMEDReader::GetFieldsTreeArrayName(int index)
405 return this->FieldSelection->GetArrayName(index);
408 //------------------------------------------------------------------------------
409 int vtkMEDReader::GetFieldsTreeArrayStatus(const char* name)
411 return this->FieldSelection->ArrayIsEnabled(name);
414 //------------------------------------------------------------------------------
415 void vtkMEDReader::SetFieldsStatus(const char* name, int status)
417 if (this->GetFieldsTreeArrayStatus(name) != status)
421 this->FieldSelection->EnableArray(name);
425 this->FieldSelection->DisableArray(name);
431 //------------------------------------------------------------------------------
432 int vtkMEDReader::GetNumberOfTimesFlagsArrays()
434 return this->TimeFlagSelection->GetNumberOfArrays();
437 //------------------------------------------------------------------------------
438 const char* vtkMEDReader::GetTimesFlagsArrayName(int index)
440 return this->TimeFlagSelection->GetArrayName(index);
443 //------------------------------------------------------------------------------
444 int vtkMEDReader::GetTimesFlagsArrayStatus(const char* name)
446 return this->TimeFlagSelection->ArrayIsEnabled(name);
449 //------------------------------------------------------------------------------
450 void vtkMEDReader::SetTimesFlagsStatus(const char* name, int status)
452 if (this->GetTimesFlagsArrayStatus(name) != status)
456 this->TimeFlagSelection->EnableArray(name);
460 this->TimeFlagSelection->DisableArray(name);
466 void vtkMEDReader::UpdateSIL(vtkInformation* /*request*/, vtkInformation * /*info*/)
470 std::string meshName(this->Internal->Tree.getActiveMeshName());
471 if(!this->Internal->SIL || meshName!=this->Internal->DftMeshName)
473 vtkMutableDirectedGraph *sil(vtkMutableDirectedGraph::New());
475 if(this->Internal->SIL)
476 this->Internal->SIL->Delete();
477 this->Internal->SIL=sil;
478 this->Internal->DftMeshName=meshName;
483 * The returned string is the name of the mesh activated which groups and families are in \a sil.
485 std::string vtkMEDReader::BuildSIL(vtkMutableDirectedGraph* sil)
488 return std::string();
490 vtkSmartPointer<vtkVariantArray> childEdge(vtkSmartPointer<vtkVariantArray>::New());
491 childEdge->InsertNextValue(0);
492 vtkSmartPointer<vtkVariantArray> crossEdge(vtkSmartPointer<vtkVariantArray>::New());
493 crossEdge->InsertNextValue(1);
494 // CrossEdge is an edge linking hierarchies.
495 vtkUnsignedCharArray* crossEdgesArray=vtkUnsignedCharArray::New();
496 crossEdgesArray->SetName("CrossEdges");
497 sil->GetEdgeData()->AddArray(crossEdgesArray);
498 crossEdgesArray->Delete();
499 std::vector<std::string> names;
500 // Now build the hierarchy.
501 vtkIdType rootId=sil->AddVertex();
502 names.push_back("SIL");
503 // Add global fields root
504 vtkIdType fieldsRoot(sil->AddChild(rootId,childEdge));
505 names.push_back("FieldsStatusTree");
506 this->Internal->Tree.feedSIL(sil,fieldsRoot,childEdge,names);
507 vtkIdType meshesFamsGrpsRoot(sil->AddChild(rootId,childEdge));
508 names.push_back("MeshesFamsGrps");
509 std::string dftMeshName(this->Internal->Tree.feedSILForFamsAndGrps(sil,meshesFamsGrpsRoot,childEdge,names));
510 // This array is used to assign names to nodes.
511 vtkStringArray *namesArray(vtkStringArray::New());
512 namesArray->SetName("Names");
513 namesArray->SetNumberOfTuples(sil->GetNumberOfVertices());
514 sil->GetVertexData()->AddArray(namesArray);
515 namesArray->Delete();
516 std::vector<std::string>::const_iterator iter;
518 for(cc=0, iter=names.begin(); iter!=names.end(); ++iter, ++cc)
519 namesArray->SetValue(cc,(*iter).c_str());
523 double vtkMEDReader::PublishTimeStepsIfNeeded(vtkInformation *outInfo, bool& isUpdated)
529 std::vector<double> tsteps;
530 if(!this->IsStdOrMode)
531 tsteps=this->Internal->Tree.getTimeSteps(lev0,this->Internal->TK);
533 { tsteps.resize(1); tsteps[0]=0.; }
535 if(lev0!=this->Internal->LastLev0)
539 timeRange[0]=tsteps.front();
540 timeRange[1]=tsteps.back();
541 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),&tsteps[0],(int)tsteps.size());
542 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),timeRange,2);
543 this->Internal->LastLev0=lev0;
545 return tsteps.front();
548 void vtkMEDReader::FillMultiBlockDataSetInstance(vtkMultiBlockDataSet *output, double reqTS, ExportedTinyInfo *internalInfo)
550 if( !this->Internal )
552 vtkDataSet *ret(RetrieveDataSetAtTime(reqTS,internalInfo));
553 output->SetBlock(0,ret);
557 vtkDataSet *vtkMEDReader::RetrieveDataSetAtTime(double reqTS, ExportedTinyInfo *internalInfo)
559 if( !this->Internal )
561 std::string meshName;
562 vtkDataSet *ret(this->Internal->Tree.buildVTKInstance(this->IsStdOrMode,reqTS,meshName,this->Internal->TK,!this->RemoveDebugArrays,internalInfo));
563 if(this->GenerateVect)
565 vtkGenerateVectors::Operate(ret->GetPointData());
566 vtkGenerateVectors::Operate(ret->GetCellData());
567 vtkGenerateVectors::Operate(ret->GetFieldData());
568 // The operations above have potentially created new arrays -> This breaks the optimization of StaticMesh that expects the same field arrays over time.
569 // To enforce the cache recomputation declare modification of mesh.
570 //vtkGenerateVectors::ChangeMeshTimeToUpdateCache(ret);
575 void vtkMEDReader::PrintSelf(ostream& os, vtkIndent indent)
577 this->Superclass::PrintSelf(os, indent);