src/Plugins/MEDReader/IO/vtkMEDReader.cxx

   1 // Copyright (C) 2010-2014  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 : Anthony Geay
  20
  21 #include "vtkMEDReader.h"
  22 #include "vtkGenerateVectors.h"
  23
  24 #include "vtkMultiBlockDataSet.h"
  25 #include "vtkInformation.h"
  26 #include "vtkDataSetAttributes.h"
  27 #include "vtkStringArray.h"
  28 #include "vtkMutableDirectedGraph.h"
  29 #include "vtkInformationStringKey.h"
  30 //
  31 #include "vtkUnsignedCharArray.h"
  32 #include "vtkInformationVector.h"
  33 #include "vtkSmartPointer.h"
  34 #include "vtkVariantArray.h"
  35 #include "vtkExecutive.h"
  36 #include "vtkStreamingDemandDrivenPipeline.h"
  37 #include "vtkMultiTimeStepAlgorithm.h"
  38 #include "vtkUnstructuredGrid.h"
  39 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
  40 #include "vtkQuadratureSchemeDefinition.h"
  41 #include "vtkPointData.h"
  42 #include "vtkCellData.h"
  43 #include "vtkCellType.h"
  44 #include "vtkCellArray.h"
  45 #include "vtkDoubleArray.h"
  46 #include "vtkObjectFactory.h"
  47
  48 #include "MEDFileFieldRepresentationTree.hxx"
  49
  50 #include <string>
  51 #include <vector>
  52 #include <sstream>
  53 #include <algorithm>
  54
  55 /*!
  56  * This class stores properties in loading state mode (pvsm) when the MED file has not been read yet.
  57  * The file is not read beacause FileName has not been informed yet ! So this class stores properties of vtkMEDReader instance that
  58  * owns it and wait the vtkMEDReader::SetFileName to apply properties afterwards.
  59  */
  60 class PropertyKeeper
  61 {
  62 public:
  63   PropertyKeeper(vtkMEDReader *master):_master(master),IsGVActivated(false),GVValue(0),IsCMActivated(false),CMValue(0) { }
  64   void assignPropertiesIfNeeded();
  65   bool arePropertiesOnTreeToSetAfter() const;
  66   //
  67   void pushFieldStatusEntry(const char* name, int status);
  68   void pushGenerateVectorsValue(int value);
  69   void pushChangeModeValue(int value);
  70   void pushTimesFlagsStatusEntry(const char* name, int status);
  71 protected:
  72   // pool of pairs to assign in SetFieldsStatus if needed. The use case is the load using pvsm.
  73   std::vector< std::pair<std::string,int> > SetFieldsStatusPairs;
  74   // generate vector
  75   bool IsGVActivated;
  76   int GVValue;
  77   // change mode
  78   bool IsCMActivated;
  79   int CMValue;
  80   //
  81   std::vector< std::pair<std::string,int> > TimesFlagsStatusPairs;
  82   vtkMEDReader *_master;
  83 };
  84
  85 void PropertyKeeper::assignPropertiesIfNeeded()
  86 {
  87   if(!this->SetFieldsStatusPairs.empty())
  88     {
  89       for(std::vector< std::pair<std::string,int> >::const_iterator it=this->SetFieldsStatusPairs.begin();it!=this->SetFieldsStatusPairs.end();it++)
  90         _master->SetFieldsStatus((*it).first.c_str(),(*it).second);
  91       this->SetFieldsStatusPairs.clear();
  92     }
  93   if(!this->TimesFlagsStatusPairs.empty())
  94     {
  95       for(std::vector< std::pair<std::string,int> >::const_iterator it=this->TimesFlagsStatusPairs.begin();it!=this->TimesFlagsStatusPairs.end();it++)
  96         _master->SetTimesFlagsStatus((*it).first.c_str(),(*it).second);
  97       this->TimesFlagsStatusPairs.clear();
  98     }
  99   if(this->IsGVActivated)
 100     {
 101       _master->GenerateVectors(this->GVValue);
 102       this->IsGVActivated=false;
 103     }
 104   if(this->IsCMActivated)
 105     {
 106       _master->ChangeMode(this->CMValue);
 107       this->IsCMActivated=false;
 108     }
 109 }
 110
 111 void PropertyKeeper::pushFieldStatusEntry(const char* name, int status)
 112 {
 113   bool found(false);
 114   for(std::vector< std::pair<std::string,int> >::const_iterator it=this->SetFieldsStatusPairs.begin();it!=this->SetFieldsStatusPairs.end() && !found;it++)
 115     found=(*it).first==name;
 116   if(!found)
 117     this->SetFieldsStatusPairs.push_back(std::pair<std::string,int>(name,status));
 118 }
 119
 120 void PropertyKeeper::pushTimesFlagsStatusEntry(const char* name, int status)
 121 {
 122   bool found(false);
 123   for(std::vector< std::pair<std::string,int> >::const_iterator it=this->TimesFlagsStatusPairs.begin();it!=this->TimesFlagsStatusPairs.end() && !found;it++)
 124     found=(*it).first==name;
 125   if(!found)
 126     this->TimesFlagsStatusPairs.push_back(std::pair<std::string,int>(name,status));
 127 }
 128
 129 void PropertyKeeper::pushGenerateVectorsValue(int value)
 130 {
 131   this->IsGVActivated=true;
 132   this->GVValue=value;
 133 }
 134
 135 void PropertyKeeper::pushChangeModeValue(int value)
 136 {
 137   this->IsCMActivated=true;
 138   this->CMValue=value;
 139 }
 140
 141 bool PropertyKeeper::arePropertiesOnTreeToSetAfter() const
 142 {
 143   return !SetFieldsStatusPairs.empty();
 144 }
 145
 146 class vtkMEDReader::vtkMEDReaderInternal
 147 {
 148
 149 public:
 150   vtkMEDReaderInternal(vtkMEDReader *master):TK(0),IsMEDOrSauv(true),IsStdOrMode(false),GenerateVect(false),SIL(0),LastLev0(-1),FirstCall0(2),PK(master)
 151   {
 152   }
 153
 154   bool PluginStart0()
 155   {
 156     if(FirstCall0==0)
 157       return false;
 158     FirstCall0--;
 159     return true;
 160   }
 161
 162   ~vtkMEDReaderInternal()
 163   {
 164     if(this->SIL)
 165       this->SIL->Delete();
 166   }
 167 public:
 168   MEDFileFieldRepresentationTree Tree;
 169   TimeKeeper TK;
 170   std::string FileName;
 171   //when true the file is MED file. when false it is a Sauv file
 172   bool IsMEDOrSauv;
 173   //when false -> std, true -> mode. By default std (false).
 174   bool IsStdOrMode;
 175   //when false -> do nothing. When true cut off or extend to nbOfCompo=3 vector arrays.
 176   bool GenerateVect;
 177   std::string DftMeshName;
 178   // Store the vtkMutableDirectedGraph that represents links between family, groups and cell types
 179   vtkMutableDirectedGraph* SIL;
 180   // store the lev0 id in Tree corresponding to the TIME_STEPS in the pipeline.
 181   int LastLev0;
 182   // The property keeper is usable only in pvsm mode.
 183   PropertyKeeper PK;
 184 private:
 185   unsigned char FirstCall0;
 186 };
 187
 188 vtkStandardNewMacro(vtkMEDReader);
 189
 190 vtkMEDReader::vtkMEDReader():Internal(new vtkMEDReaderInternal(this))
 191 {
 192   this->SetNumberOfInputPorts(0);
 193   this->SetNumberOfOutputPorts(1);
 194 }
 195
 196 vtkMEDReader::~vtkMEDReader()
 197 {
 198   delete this->Internal;
 199 }
 200
 201 void vtkMEDReader::Reload(int a)
 202 {
 203   if(a==0)
 204     return;
 205   std::cerr << "vtkMEDReader::Reload" << a << std::endl;
 206   std::string fName((const char *)this->GetFileName());
 207   delete this->Internal;
 208   this->Internal=new vtkMEDReaderInternal(this);
 209   this->SetFileName(fName.c_str());
 210 }
 211
 212 void vtkMEDReader::GenerateVectors(int val)
 213 {
 214   if(this->Internal->FileName.empty())
 215     {//pvsm mode
 216       this->Internal->PK.pushGenerateVectorsValue(val);
 217       return ;
 218     }
 219   //not pvsm mode (general case)
 220   bool val2((bool)val);
 221   if(val2!=this->Internal->GenerateVect)
 222     {
 223       this->Internal->GenerateVect=val2;
 224       this->Modified();
 225     }
 226 }
 227
 228 void vtkMEDReader::ChangeMode(int newMode)
 229 {
 230   if(this->Internal->FileName.empty())
 231     {//pvsm mode
 232       this->Internal->PK.pushChangeModeValue(newMode);
 233       return ;
 234     }
 235   //not pvsm mode (general case)
 236   this->Internal->IsStdOrMode=newMode!=0;
 237   //std::cerr << "vtkMEDReader::ChangeMode : " << this->Internal->IsStdOrMode << std::endl;
 238   this->Modified();
 239 }
 240
 241 const char *vtkMEDReader::GetSeparator()
 242 {
 243   return MEDFileFieldRepresentationLeavesArrays::ZE_SEP;
 244 }
 245
 246 void vtkMEDReader::SetFileName(const char *fname)
 247 {
 248   try
 249     {
 250       this->Internal->FileName=fname;
 251       std::size_t pos(this->Internal->FileName.find_last_of('.'));
 252       if(pos!=std::string::npos)
 253         {
 254           std::string ext(this->Internal->FileName.substr(pos));
 255           if(ext.find("sauv")!=std::string::npos)
 256             this->Internal->IsMEDOrSauv=false;
 257         }
 258       if(this->Internal->Tree.getNumberOfLeavesArrays()==0)
 259         {
 260           int iPart(-1),nbOfParts(-1);
 261 #ifdef MEDREADER_USE_MPI
 262           MPI_Comm_rank(MPI_COMM_WORLD,&iPart);
 263           MPI_Comm_size(MPI_COMM_WORLD,&nbOfParts);
 264 #endif
 265           this->Internal->Tree.loadMainStructureOfFile(this->Internal->FileName.c_str(),this->Internal->IsMEDOrSauv,iPart,nbOfParts);
 266           if(!this->Internal->PK.arePropertiesOnTreeToSetAfter())
 267             this->Internal->Tree.activateTheFirst();//This line manually initialize the status of server (this) with the remote client.
 268           this->Internal->TK.setMaxNumberOfTimeSteps(this->Internal->Tree.getMaxNumberOfTimeSteps());
 269         }
 270       this->Modified();
 271       this->Internal->PK.assignPropertiesIfNeeded();
 272     }
 273   catch(INTERP_KERNEL::Exception& e)
 274     {
 275       delete this->Internal;
 276       this->Internal=0;
 277       std::ostringstream oss;
 278       oss << "Exception has been thrown in vtkMEDReader::SetFileName : " << e.what() << std::endl;
 279       if(this->HasObserver("ErrorEvent") )
 280         this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
 281       else
 282         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
 283       vtkObject::BreakOnError();
 284     }
 285 }
 286
 287 char *vtkMEDReader::GetFileName()
 288 {
 289   return const_cast<char *>(this->Internal->FileName.c_str());
 290 }
 291
 292 int vtkMEDReader::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 293 {
 294   //std::cerr << "########################################## RequestInformation ##########################################" << std::endl;
 295   if(!this->Internal)
 296     return 0;
 297   try
 298     {
 299       vtkInformation *outInfo(outputVector->GetInformationObject(0));
 300       outInfo->Set(vtkStreamingDemandDrivenPipeline::MAXIMUM_NUMBER_OF_PIECES(),-1);
 301       outInfo->Set(vtkDataObject::DATA_TYPE_NAME(),"vtkMultiBlockDataSet");
 302       this->UpdateSIL(outInfo);
 303       //
 304       bool dummy(false);
 305       this->PublishTimeStepsIfNeeded(outInfo,dummy);
 306     }
 307   catch(INTERP_KERNEL::Exception& e)
 308     {
 309       std::cerr << "Exception has been thrown in vtkMEDReader::RequestInformation : " << e.what() << std::endl;
 310       return 0;
 311     }
 312   return 1;
 313 }
 314
 315 int vtkMEDReader::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 316 {
 317   //std::cerr << "########################################## RequestData        ##########################################";
 318   if(!this->Internal)
 319     return 0;
 320   try
 321     {
 322       vtkInformation *outInfo(outputVector->GetInformationObject(0));
 323       vtkMultiBlockDataSet *output(vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
 324       bool isUpdated(false);
 325       double reqTS(0.);
 326       if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
 327         reqTS=outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
 328       //std::cerr << reqTS << std::endl;
 329       this->FillMultiBlockDataSetInstance(output,reqTS);
 330       output->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(),reqTS);
 331       this->UpdateSIL(outInfo);
 332       //this->UpdateProgress((float) progress/((float) maxprogress-1));
 333     }
 334   catch(INTERP_KERNEL::Exception& e)
 335     {
 336       std::cerr << "Exception has been thrown in vtkMEDReader::RequestInformation : " << e.what() << std::endl;
 337       return 0;
 338     }
 339   return 1;
 340 }
 341
 342 void vtkMEDReader::SetFieldsStatus(const char* name, int status)
 343 {
 344   //std::cerr << "vtkMEDReader::SetFieldsStatus(" << name << "," << status << ") called !" << std::endl;
 345   if(this->Internal->FileName.empty())
 346     {//pvsm mode
 347       this->Internal->PK.pushFieldStatusEntry(name,status);
 348       return ;
 349     }
 350   //not pvsm mode (general case)
 351   try
 352     {
 353       this->Internal->Tree.changeStatusOfAndUpdateToHaveCoherentVTKDataSet(this->Internal->Tree.getIdHavingZeName(name),status);
 354       if(std::string(name)==GetFieldsTreeArrayName(GetNumberOfFieldsTreeArrays()-1))
 355         if(!this->Internal->PluginStart0())
 356           this->Modified();
 357     }
 358   catch(INTERP_KERNEL::Exception& e)
 359     {
 360       if(!this->Internal->FileName.empty())
 361         {
 362           std::cerr << "vtkMEDReader::SetFieldsStatus error : " << e.what() << " *** WITH STATUS=" << status << endl;
 363           return ;
 364         }
 365     }
 366 }
 367
 368 int vtkMEDReader::GetNumberOfFieldsTreeArrays()
 369 {
 370   if(!this->Internal)
 371     return 0;
 372   int ret(this->Internal->Tree.getNumberOfLeavesArrays());
 373   //std::cerr << "vtkMEDReader::GetNumberOfFieldsTreeArrays called ! " << ret << std::endl;
 374   return ret;
 375 }
 376
 377 const char *vtkMEDReader::GetFieldsTreeArrayName(int index)
 378 {
 379   return this->Internal->Tree.getNameOfC(index);
 380   //std::cerr << "vtkMEDReader::GetFieldsTreeArrayName(" << index << ") called ! " << ret << std::endl;
 381 }
 382
 383 int vtkMEDReader::GetFieldsTreeArrayStatus(const char *name)
 384 {
 385   int zeId(this->Internal->Tree.getIdHavingZeName(name));
 386   int ret(this->Internal->Tree.getStatusOf(zeId));
 387   return ret;
 388 }
 389
 390 void vtkMEDReader::SetTimesFlagsStatus(const char *name, int status)
 391 {
 392   if(this->Internal->FileName.empty())
 393     {//pvsm mode
 394       this->Internal->PK.pushTimesFlagsStatusEntry(name,status);
 395       return ;
 396     }
 397   //not pvsm mode (general case)
 398   int pos(0);
 399   std::istringstream iss(name); iss >> pos;
 400   this->Internal->TK.getTimesFlagArray()[pos].first=(bool)status;
 401   if(pos==this->Internal->TK.getTimesFlagArray().size()-1)
 402     if(!this->Internal->PluginStart0())
 403       {
 404         this->Modified();
 405         //this->Internal->TK.printSelf(std::cerr);
 406       }
 407 }
 408
 409 int vtkMEDReader::GetNumberOfTimesFlagsArrays()
 410 {
 411   if(!this->Internal)
 412     return 0;
 413   return (int)this->Internal->TK.getTimesFlagArray().size();
 414 }
 415
 416 const char *vtkMEDReader::GetTimesFlagsArrayName(int index)
 417 {
 418   return this->Internal->TK.getTimesFlagArray()[index].second.c_str();
 419 }
 420
 421 int vtkMEDReader::GetTimesFlagsArrayStatus(const char *name)
 422 {
 423   int pos(0);
 424   std::istringstream iss(name); iss >> pos;
 425   return (int)this->Internal->TK.getTimesFlagArray()[pos].first;
 426 }
 427
 428 void vtkMEDReader::UpdateSIL(vtkInformation *info)
 429 {
 430   if(!this->Internal)
 431       return ;
 432   vtkMutableDirectedGraph *sil(vtkMutableDirectedGraph::New());
 433   std::string meshName(this->BuildSIL(sil));
 434   if(meshName!=this->Internal->DftMeshName)
 435     {
 436       if(this->Internal->SIL)
 437         this->Internal->SIL->Delete();
 438       this->Internal->SIL=sil;
 439       this->Internal->DftMeshName=meshName;
 440       info->Set(vtkDataObject::SIL(),this->Internal->SIL);
 441       //request->AppendUnique(vtkExecutive::KEYS_TO_COPY(),vtkDataObject::SIL());
 442     }
 443   else
 444     {
 445       sil->Delete();
 446     }
 447 }
 448
 449 /*!
 450  * The returned string is the name of the mesh activated which groups and families are in \a sil.
 451  */
 452 std::string vtkMEDReader::BuildSIL(vtkMutableDirectedGraph* sil)
 453 {
 454   sil->Initialize();
 455   vtkSmartPointer<vtkVariantArray> childEdge(vtkSmartPointer<vtkVariantArray>::New());
 456   childEdge->InsertNextValue(0);
 457   vtkSmartPointer<vtkVariantArray> crossEdge(vtkSmartPointer<vtkVariantArray>::New());
 458   crossEdge->InsertNextValue(1);
 459   // CrossEdge is an edge linking hierarchies.
 460   vtkUnsignedCharArray* crossEdgesArray=vtkUnsignedCharArray::New();
 461   crossEdgesArray->SetName("CrossEdges");
 462   sil->GetEdgeData()->AddArray(crossEdgesArray);
 463   crossEdgesArray->Delete();
 464   std::vector<std::string> names;
 465   // Now build the hierarchy.
 466   vtkIdType rootId=sil->AddVertex();
 467   names.push_back("SIL");
 468   // Add global fields root
 469   vtkIdType fieldsRoot(sil->AddChild(rootId,childEdge));
 470   names.push_back("FieldsStatusTree");
 471   this->Internal->Tree.feedSIL(sil,fieldsRoot,childEdge,names);
 472   vtkIdType meshesFamsGrpsRoot(sil->AddChild(rootId,childEdge));
 473   names.push_back("MeshesFamsGrps");
 474   std::string dftMeshName(this->Internal->Tree.feedSILForFamsAndGrps(sil,meshesFamsGrpsRoot,childEdge,names));
 475   // This array is used to assign names to nodes.
 476   vtkStringArray *namesArray(vtkStringArray::New());
 477   namesArray->SetName("Names");
 478   namesArray->SetNumberOfTuples(sil->GetNumberOfVertices());
 479   sil->GetVertexData()->AddArray(namesArray);
 480   namesArray->Delete();
 481   std::vector<std::string>::const_iterator iter;
 482   vtkIdType cc;
 483   for(cc=0, iter=names.begin(); iter!=names.end(); ++iter, ++cc)
 484     namesArray->SetValue(cc,(*iter).c_str());
 485   return dftMeshName;
 486 }
 487
 488 double vtkMEDReader::PublishTimeStepsIfNeeded(vtkInformation *outInfo, bool& isUpdated)
 489 {
 490   int lev0(-1);
 491   std::vector<double> tsteps;
 492   if(!this->Internal->IsStdOrMode)
 493     tsteps=this->Internal->Tree.getTimeSteps(lev0,this->Internal->TK);
 494   else
 495     { tsteps.resize(1); tsteps[0]=0.; }
 496   isUpdated=false;
 497   if(lev0!=this->Internal->LastLev0)
 498     {
 499       isUpdated=true;
 500       double timeRange[2];
 501       timeRange[0]=tsteps.front();
 502       timeRange[1]=tsteps.back();
 503       outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),&tsteps[0],tsteps.size());
 504       outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),timeRange,2);
 505       this->Internal->LastLev0=lev0;
 506     }
 507   return tsteps.front();
 508 }
 509
 510 void vtkMEDReader::FillMultiBlockDataSetInstance(vtkMultiBlockDataSet *output, double reqTS)
 511 {
 512   std::string meshName;
 513   vtkDataSet *ret(this->Internal->Tree.buildVTKInstance(this->Internal->IsStdOrMode,reqTS,meshName,this->Internal->TK));
 514   if(this->Internal->GenerateVect)
 515     {
 516       vtkGenerateVectors::Operate(ret->GetPointData());
 517       vtkGenerateVectors::Operate(ret->GetCellData());
 518       vtkGenerateVectors::Operate(ret->GetFieldData());
 519     }
 520   output->SetBlock(0,ret);
 521   ret->Delete();
 522 }
 523
 524 void vtkMEDReader::PrintSelf(ostream& os, vtkIndent indent)
 525 {
 526   this->Superclass::PrintSelf(os, indent);
 527 }