src/Plugins/MedReader/IO/vtkMedReader.cxx

   1 // Copyright (C) 2010-2013  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.
   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
  20 #include "vtkMedReader.h"
  21
  22 #include "vtkMedFile.h"
  23 #include "vtkMedDriver.h"
  24 #include "vtkMedUtilities.h"
  25 #include "vtkMedFamily.h"
  26 #include "vtkMedGroup.h"
  27 #include "vtkMedMesh.h"
  28 #include "vtkMedUnstructuredGrid.h"
  29 #include "vtkMedCurvilinearGrid.h"
  30 #include "vtkMedRegularGrid.h"
  31 #include "vtkMedEntityArray.h"
  32 #include "vtkMedField.h"
  33 #include "vtkMedFieldStep.h"
  34 #include "vtkMedFieldOverEntity.h"
  35 #include "vtkMedProfile.h"
  36 #include "vtkMedIntArray.h"
  37 #include "vtkMedLocalization.h"
  38 #include "vtkMedFamilyOnEntity.h"
  39 #include "vtkMedFieldOnProfile.h"
  40 #include "vtkMedSelection.h"
  41 #include "vtkMedFamilyOnEntityOnProfile.h"
  42 #include "vtkMedLink.h"
  43 #include "vtkMedInterpolation.h"
  44 #include "vtkMedStructElement.h"
  45 #include "vtkMedConstantAttribute.h"
  46
  47 #include "vtkObjectFactory.h"
  48 #include "vtkMutableDirectedGraph.h"
  49 #include "vtkStringArray.h"
  50 #include "vtkDataSetAttributes.h"
  51 #include "vtkUnsignedCharArray.h"
  52 #include "vtkInformation.h"
  53 #include "vtkInformationVector.h"
  54 #include "vtkSmartPointer.h"
  55 #include "vtkVariantArray.h"
  56 #include "vtkMultiBlockDataSet.h"
  57 #include "vtkExecutive.h"
  58 #include "vtkStreamingDemandDrivenPipeline.h"
  59 #include "vtkMultiTimeStepAlgorithm.h"
  60 #include "vtkUnstructuredGrid.h"
  61 #include "vtkMath.h"
  62 #include "vtkPointData.h"
  63 #include "vtkCellData.h"
  64 #include "vtkFieldData.h"
  65 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
  66 #include "vtkQuadratureSchemeDefinition.h"
  67 #include "vtkCellType.h"
  68 #include "vtkCellArray.h"
  69 #include "vtkDoubleArray.h"
  70 #include "vtkConfigure.h"
  71 #include "vtkMultiProcessController.h"
  72 #include "vtkCommunicator.h"
  73
  74 #include "vtkSMDoubleVectorProperty.h"
  75 #include "vtkInformationDataObjectKey.h"
  76
  77 #include <deque>
  78 #include <map>
  79 #include <string>
  80 #include <sstream>
  81 #include <vector>
  82 #include <list>
  83 #include <set>
  84 #include <algorithm>
  85 using namespace std;
  86
  87 struct VTKField
  88 {
  89   vtkSmartPointer<vtkDataArray> DataArray;
  90   vtkSmartPointer<vtkDataArray> Vectors;
  91   vtkSmartPointer<vtkIdTypeArray> QuadratureIndexArray;
  92 };
  93
  94 class vtkMedListOfFieldSteps : public std::list<vtkMedFieldStep*>
  95 {};
  96
  97 typedef int LocalizationKey;
  98
  99 class vtkMedReader::vtkMedReaderInternal
 100 {
 101 public:
 102   int NumberOfPieces;
 103   int CurrentPieceNumber;
 104   int GhostLevel;
 105   double UpdateTimeStep;
 106   vtkTimeStamp FileNameMTime;
 107   vtkTimeStamp MetaDataMTime;
 108   vtkTimeStamp GroupSelectionMTime;
 109   vtkTimeStamp FamilySelectionMTime;
 110   int SILUpdateStamp;
 111   int RealAnimationMode;
 112   vtkMedSelection* Families;
 113
 114   // this stores the aggregation of all compute steps from
 115   // both meshes and fields.
 116   std::map<med_float, std::set<med_int> > GlobalComputeStep;
 117
 118   // Store the vtkMutableDirectedGraph that represents links between family, groups and cell types
 119   vtkMutableDirectedGraph* SIL;
 120
 121   // this map is used to keep clean data sets in the cache, without any field.
 122   // for each support, store the vtkDataSet
 123   map<vtkMedFamilyOnEntityOnProfile*, vtkSmartPointer<vtkDataSet> > DataSetCache;
 124
 125   // this is the current dataset for the given support.
 126   map<vtkMedFamilyOnEntityOnProfile*, vtkDataSet*> CurrentDataSet;
 127
 128   // for each support, cache the VTK arrays that correspond to a given field at the given step.
 129   map<vtkMedFamilyOnEntityOnProfile*, map<vtkMedFieldOnProfile*, VTKField> > FieldCache;
 130   //map<vtkMedFamilyOnEntity*, map<vtkMedFieldOnProfile*, bool> > FieldMatchCache;
 131
 132   // This list keep tracks of all the currently selected supports
 133   set<vtkMedFamilyOnEntityOnProfile*> UsedSupports;
 134
 135   // this map keeps for each support, the quadrature offset array so that
 136   // different fields on the same support can use
 137   // the same offset array, provided they use the same gauss points definitions
 138   map<vtkMedFamilyOnEntityOnProfile*,
 139       map<LocalizationKey, vtkSmartPointer<vtkIdTypeArray> > >
 140       QuadratureOffsetCache;
 141
 142   map<vtkMedFamilyOnEntityOnProfile*,
 143       map<vtkMedFieldOnProfile*, LocalizationKey> > QuadOffsetKey;
 144
 145   std::map<std::string, vtkSmartPointer<vtkMedFile> > MedFiles;
 146
 147   vtkMedReaderInternal()
 148   {
 149     this->SIL=vtkMutableDirectedGraph::New();
 150     this->SILUpdateStamp=-1;
 151     this->RealAnimationMode=vtkMedReader::PhysicalTime;
 152     this->Families=vtkMedSelection::New();
 153     this->FamilySelectionMTime.Modified();
 154     this->GroupSelectionMTime.Modified();
 155   }
 156   ~vtkMedReaderInternal()
 157   {
 158     this->SIL->Delete();
 159     this->Families->Delete();
 160   }
 161
 162   void ClearSupportCache(vtkMedFamilyOnEntityOnProfile* foep)
 163   {
 164     //this->Med2VTKPointIndex.erase(foep);
 165     this->QuadratureOffsetCache.erase(foep);
 166     //this->FieldMatchCache.erase(foe);
 167     this->FieldCache.erase(foep);
 168     this->CurrentDataSet.erase(foep);
 169     this->DataSetCache.erase(foep);
 170   }
 171
 172   vtkIdType GetChild(vtkIdType parent, const vtkStdString childName)
 173   {
 174     vtkStringArray* names=vtkStringArray::SafeDownCast(
 175         this->SIL->GetVertexData()->GetArray("Names"));
 176     if(names==NULL)
 177       return -1;
 178     vtkIdType nedges=this->SIL->GetOutDegree(parent);
 179     for(vtkIdType id=0; id<nedges; id++)
 180       {
 181       vtkOutEdgeType edge=this->SIL->GetOutEdge(parent, id);
 182       if(names->GetValue(edge.Target)==childName)
 183         return edge.Target;
 184       }
 185     return -1;
 186   }
 187
 188   vtkIdType GetGroupId(const char* key)
 189   {
 190     std::string meshname, celltypename, groupname;
 191     vtkMedUtilities::SplitGroupKey(key, meshname, celltypename, groupname);
 192     vtkIdType root=GetChild(0, "SIL");
 193     if(root==-1)
 194       return -1;
 195     vtkIdType mesh=GetChild(root, meshname);
 196     if(mesh==-1)
 197       return -1;
 198     vtkIdType type=GetChild(mesh, celltypename);
 199     if(type==-1)
 200       return -1;
 201     return GetChild(type, groupname);
 202
 203   }
 204
 205 };
 206
 207 //vtkCxxRevisionMacro(vtkMedReader, "$Revision$");
 208 vtkStandardNewMacro(vtkMedReader);
 209
 210 vtkMedReader::vtkMedReader()
 211 {
 212   this->FileName=NULL;
 213   this->SetNumberOfInputPorts(0);
 214   this->PointFields=vtkMedSelection::New();
 215   this->CellFields=vtkMedSelection::New();
 216   this->QuadratureFields=vtkMedSelection::New();
 217   this->ElnoFields=vtkMedSelection::New();
 218   this->Entities=vtkMedSelection::New();
 219   this->Groups=vtkMedSelection::New();
 220   this->Frequencies=vtkMedSelection::New();
 221   this->AnimationMode=Default;
 222   this->TimeIndexForIterations=0;
 223   this->CacheStrategy=CacheGeometry;
 224   this->Internal=new vtkMedReaderInternal;
 225   this->TimePrecision=0.00001;
 226   this->AvailableTimes=vtkDoubleArray::New();
 227   this->GenerateVectors = 0;
 228 }
 229
 230 vtkMedReader::~vtkMedReader()
 231 {
 232   this->SetFileName(NULL);
 233   this->PointFields->Delete();
 234   this->CellFields->Delete();
 235   this->QuadratureFields->Delete();
 236   this->ElnoFields->Delete();
 237   this->Entities->Delete();
 238   this->Groups->Delete();
 239   this->Frequencies->Delete();
 240   delete this->Internal;
 241   this->AvailableTimes->Delete();
 242 }
 243
 244 int vtkMedReader::GetSILUpdateStamp()
 245 {
 246   return this->Internal->SILUpdateStamp;
 247 }
 248
 249 void vtkMedReader::SetFileName(const char* fname)
 250 {
 251   int modified=0;
 252   if(fname==this->FileName)
 253     return;
 254   if(fname&&this->FileName&&!strcmp(fname, this->FileName))
 255     return;
 256   modified=1;
 257   if(this->FileName)
 258     delete[] this->FileName;
 259   if (fname)
 260     {
 261     size_t fnl=strlen(fname)+1;
 262     char* dst=new char[fnl];
 263     const char* src=fname;
 264     this->FileName=dst;
 265     do
 266       {
 267       *dst++=*src++;
 268       }
 269     while (--fnl);
 270     }
 271   else
 272     {
 273     this->FileName=0;
 274     }
 275   if (modified)
 276     {
 277     this->Modified();
 278     this->Internal->MedFiles.clear();
 279     this->Internal->FileNameMTime.Modified();
 280     }
 281 }
 282
 283 int vtkMedReader::CanReadFile(const char* fname)
 284 {
 285   // the factory give a driver only when it can read the file version,
 286   // or it returns a NULL pointer.
 287   vtkSmartPointer<vtkMedFile> file=vtkSmartPointer<vtkMedFile>::New();
 288   file->SetFileName(fname);
 289
 290   if(!file->CreateDriver())
 291     return 0;
 292
 293   return 1;
 294 }
 295
 296 int vtkMedReader::RequestInformation(vtkInformation *request,
 297     vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 298 {
 299   vtkInformation* outInfo = outputVector->GetInformationObject(0);
 300   outInfo->Set(vtkStreamingDemandDrivenPipeline::MAXIMUM_NUMBER_OF_PIECES(),-1);
 301
 302   if(this->Internal->MetaDataMTime <= this->Internal->FileNameMTime)
 303     {
 304     this->ClearCaches(Initialize);
 305
 306     vtkMedFile* file=vtkMedFile::New();
 307     file->SetFileName(this->FileName);
 308     this->Internal->MedFiles[this->FileName] = file;
 309     file->Delete();
 310
 311     std::list<vtkMedFile*> file_stack;
 312     file_stack.push_back(file);
 313
 314     // if the file cannot create a driver, this means that the filename is not
 315     // valid, or that I do not know how to read this file.
 316     while (file_stack.size() > 0)
 317       {
 318       vtkMedFile* file = file_stack.front();
 319       file_stack.pop_front();
 320       // This reads information from disk
 321       file->ReadInformation();
 322
 323       // add all files linked to in the current file to the files to read.
 324       for (int linkid=0; linkid<file->GetNumberOfLink(); linkid++)
 325         {
 326         vtkMedLink* link = file->GetLink(linkid);
 327         const char* filename = link->GetFullLink(file->GetFileName());
 328         if(this->Internal->MedFiles.find(filename) == this->Internal->MedFiles.end())
 329           {
 330           vtkMedFile* newfile = vtkMedFile::New();
 331           newfile->SetFileName(filename);
 332           this->Internal->MedFiles[filename] = newfile;
 333           file_stack.push_back(newfile);
 334           newfile->Delete();
 335           }
 336         }
 337       }
 338
 339     // This computes some meta information, like which field use which
 340     // support, but do not read large data from disk.
 341     this->LinkMedInfo();
 342
 343     // This computes the initial global id of each cell type.
 344     this->InitializeCellGlobalIds();
 345
 346     this->ClearSelections();
 347
 348     this->BuildSIL(this->Internal->SIL);
 349     this->Internal->SILUpdateStamp++;
 350
 351     this->GatherComputeSteps();
 352
 353     this->Internal->MetaDataMTime.Modified();
 354     }
 355
 356   outInfo->Set(vtkDataObject::SIL(), this->Internal->SIL);
 357   request->AppendUnique(vtkExecutive::KEYS_TO_COPY(),
 358                         vtkDataObject::SIL());
 359   request->AppendUnique(vtkExecutive::KEYS_TO_COPY(),
 360                         vtkMedUtilities::BLOCK_NAME());
 361   this->AdvertiseTime(outInfo);
 362   return 1;
 363 }
 364
 365 int vtkMedReader::RequestData(vtkInformation *request,
 366     vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 367 {
 368   if(this->FileName==NULL)
 369     {
 370     vtkWarningMacro( << "FileName must be set and meta data loaded");
 371     return 0;
 372     }
 373
 374   vtkInformation *info=outputVector->GetInformationObject(0);
 375
 376   vtkMultiBlockDataSet *output=vtkMultiBlockDataSet::SafeDownCast(info->Get(
 377       vtkDataObject::DATA_OBJECT()));
 378
 379   if (info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()))
 380     {
 381     this->Internal->NumberOfPieces=info->Get(
 382         vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
 383     }
 384   else
 385     {
 386     vtkMultiProcessController* controller =
 387                 vtkMultiProcessController::GetGlobalController();
 388     if(controller)
 389       {
 390       this->Internal->NumberOfPieces=controller->GetNumberOfProcesses();
 391       }
 392     else
 393       {
 394       this->Internal->NumberOfPieces = 1;
 395       }
 396     }
 397   if (info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()))
 398     {
 399     this->Internal->CurrentPieceNumber=info->Get(
 400         vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
 401     }
 402   else
 403     {
 404     this->Internal->CurrentPieceNumber=0;
 405     vtkMultiProcessController* controller =
 406                 vtkMultiProcessController::GetGlobalController();
 407     if(controller)
 408       {
 409       this->Internal->CurrentPieceNumber= controller->GetLocalProcessId();
 410       }
 411     else
 412       {
 413       this->Internal->CurrentPieceNumber=0;
 414       }
 415     }
 416
 417   if (info->Has(
 418       vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()))
 419     {
 420     this->Internal->GhostLevel=info->Get(
 421         vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
 422     }
 423   else
 424     {
 425     this->Internal->GhostLevel=0;
 426     }
 427
 428   //vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()
 429   /* if (info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()))
 430     {
 431     this->Internal->UpdateTimeStep=info->Get(
 432         vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())[0];
 433     }
 434   else
 435     {
 436     this->Internal->UpdateTimeStep=0;
 437     } */
 438
 439   if (info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
 440     {
 441     this->Internal->UpdateTimeStep=info->Get(
 442         vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
 443     }
 444   else
 445     {
 446     this->Internal->UpdateTimeStep=0;
 447     }
 448
 449   this->InitializeParallelRead();
 450   output->Initialize();
 451
 452   this->ChooseRealAnimationMode();
 453
 454   std::list<vtkMedDriver::FileOpen> openlist;
 455   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 456       fileit = this->Internal->MedFiles.begin();
 457   while(fileit != this->Internal->MedFiles.end())
 458     {
 459     vtkMedFile* file = fileit->second;
 460     openlist.push_back(vtkMedDriver::FileOpen(file->GetMedDriver()));
 461     fileit++;
 462     }
 463
 464   // clear the dataset cache of unneeded geometry
 465   this->ClearCaches(StartRequest);
 466
 467   // This call create the vtkMedSupports, but do not create the corresponding vtkDataSet;
 468   this->CreateMedSupports();
 469   this->ClearCaches(AfterCreateMedSupports);
 470   // This call creates the actual vtkDataSet that corresponds to each support
 471   int supportId = 0;
 472   int progress=0;
 473   int maxprogress=2*this->Internal->UsedSupports.size();
 474   supportId = 0;
 475   int it_counter = 0;
 476   for(set<vtkMedFamilyOnEntityOnProfile*>::iterator it=
 477       this->Internal->UsedSupports.begin(); it
 478       !=this->Internal->UsedSupports.end(); it++)
 479     {
 480     ostringstream sstr;
 481     vtkMedFamilyOnEntityOnProfile* foep = *it;
 482     sstr<<"Support : "<<vtkMedUtilities::SimplifyName(
 483         foep->GetFamilyOnEntity()->GetFamily()->GetName());
 484     this->SetProgressText(sstr.str().c_str());
 485     int doBuildSupportField = 1;
 486     it_counter++;
 487     this->BuildVTKSupport(foep, doBuildSupportField);
 488     this->UpdateProgress((float) progress/((float) maxprogress-1));
 489     progress++;
 490     supportId++;
 491     }
 492
 493   this->ClearCaches(EndBuildVTKSupports);
 494   // This call maps the fields to the supports
 495   for(set<vtkMedFamilyOnEntityOnProfile*>::iterator it=
 496       this->Internal->UsedSupports.begin(); it
 497       !=this->Internal->UsedSupports.end(); it++)
 498     {
 499     vtkMedFamilyOnEntityOnProfile* foep = *it;
 500     if((foep->GetValid() == 0) && (this->Internal->NumberOfPieces == 1))
 501       continue;
 502     ostringstream sstr;
 503     sstr<<"Loading fields on "<<vtkMedUtilities::SimplifyName(
 504         foep->GetFamilyOnEntity()->GetFamily()->GetName());
 505     this->SetProgressText(sstr.str().c_str());
 506     int doMapField = 1;
 507     this->MapFieldsOnSupport(*it, doMapField);
 508     this->UpdateProgress((float) progress/((float) maxprogress-1));
 509     progress++;
 510     supportId++;
 511     }
 512
 513   // This call clean up caches (what is actually done depends of the CacheStrategy)
 514   this->ClearCaches(EndRequest);
 515   return 1;
 516 }
 517
 518 void vtkMedReader::InitializeCellGlobalIds()
 519 {
 520   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 521       fileit = this->Internal->MedFiles.begin();
 522   while(fileit != this->Internal->MedFiles.end())
 523     {
 524     vtkMedFile* file = fileit->second;
 525     fileit++;
 526     for(int m=0; m<file->GetNumberOfMesh(); m++)
 527       {
 528       vtkMedMesh* mesh = file->GetMesh(m);
 529       med_int nstep = mesh->GetNumberOfGridStep();
 530       for(med_int stepid=0; stepid<nstep; stepid++)
 531         {
 532         vtkMedGrid* grid = mesh->GetGridStep(stepid);
 533         grid->InitializeCellGlobalIds();
 534         }
 535       }
 536     }
 537 }
 538
 539 // Method to create the filters for the MED parallel read functions
 540 // It is defined here as we have all information for initialization
 541 void vtkMedReader::InitializeParallelRead()
 542 {
 543   // If there is only one process for reading no need to enter here
 544   if (this->Internal->NumberOfPieces <= 1)
 545     {
 546     return;
 547     }
 548
 549   // FIRST: Generate filters for the cells
 550
 551   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 552   meshfit = this->Internal->MedFiles.begin();
 553   while(meshfit != this->Internal->MedFiles.end())
 554     {
 555     vtkMedFile* meshfile = meshfit->second;
 556     meshfit++;
 557     med_idt pFileID = meshfile->GetMedDriver()->GetParallelFileId();
 558
 559     for(int mid=0; mid<meshfile->GetNumberOfMesh(); mid++)
 560       {
 561       vtkMedMesh* mesh = meshfile->GetMesh(mid);
 562       for(int gid=0; gid<mesh->GetNumberOfGridStep(); gid++)
 563         {
 564         vtkMedGrid* grid = mesh->GetGridStep(gid);
 565         // read point family data and create EntityArrays
 566
 567         for(int eid=0; eid < grid->GetNumberOfEntityArray(); eid++)
 568          {
 569           vtkMedEntityArray* array = grid->GetEntityArray(eid);
 570
 571           // Next continue is to avoid to create filters for the
 572           // points, at the moment we charge the points in all nodes
 573           if (array->GetEntity().GeometryType == MED_POINT1) // !MED_NODE
 574             continue;
 575
 576           med_int nbofconstituentpervalue = vtkMedUtilities::GetNumberOfNodes(
 577                                             array->GetEntity().GeometryType);
 578           if (nbofconstituentpervalue == -1)
 579             vtkErrorMacro("Still not implemented for MED_POLYGON and MED_POLYHEDRON"); // à gerer
 580
 581           // Calculating block sizes
 582           int nEntity = array->GetNumberOfEntity();
 583           int block_size = ( nEntity / this->Internal->NumberOfPieces );
 584           med_size    start  = block_size * this->Internal->CurrentPieceNumber + 1;
 585           med_size    stride = block_size;
 586           med_size    count  = 1;
 587           med_size    blocksize = block_size;
 588           med_size    lastblocksize = (nEntity % this->Internal->NumberOfPieces);
 589           if ((this->Internal->NumberOfPieces ==
 590               this->Internal->CurrentPieceNumber+1) && (lastblocksize != 0))
 591             {
 592             blocksize += lastblocksize;
 593             stride    += lastblocksize;
 594             }
 595           lastblocksize = 0;
 596
 597           vtkMedFilter *filter = vtkMedFilter::New();
 598           filter->SetFilterSizes( start, stride, count, blocksize, lastblocksize );
 599           array->SetFilter(filter);
 600          }//entity array
 601         }// grid step
 602       }//mesh
 603     }//mesh file
 604
 605   // SECOND: Filters for the Fields
 606
 607   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 608       fieldfileit;
 609   // link the FieldOnProfile with the profiles
 610   fieldfileit = this->Internal->MedFiles.begin();
 611   while(fieldfileit != this->Internal->MedFiles.end())
 612     {
 613     vtkMedFile* fieldfile = fieldfileit->second;
 614     fieldfileit++;
 615     med_idt pFileID = fieldfile->GetMedDriver()->GetParallelFileId();
 616
 617     for(int fid=0; fid<fieldfile->GetNumberOfField(); fid++)
 618       {
 619       vtkMedField* field = fieldfile->GetField(fid);
 620
 621       if (field->GetFieldType() == vtkMedField::CellField)
 622       {
 623       for(int sid = 0; sid< field->GetNumberOfFieldStep(); sid++)
 624         {
 625         vtkMedFieldStep* step = field->GetFieldStep(sid);
 626
 627         for(int foeid = 0; foeid < step->GetNumberOfFieldOverEntity(); foeid++)
 628         // TODO : seul le premier pas de temps est dispo au debut
 629           {
 630           vtkMedFieldOverEntity* fieldOverEntity = step->GetFieldOverEntity(foeid);
 631
 632           for(int fopid = 0; fopid < fieldOverEntity->GetNumberOfFieldOnProfile(); fopid++)
 633             {
 634             vtkMedFieldOnProfile* fop = fieldOverEntity->GetFieldOnProfile(fopid);
 635             // Here implement the filters as before:
 636             // 1- Modify vtkMedFieldOnProfile to contain a filter
 637             // 2- Create the filters here only if they are on CELLs (use GetFieldType)
 638             med_int nbofconstituentpervalue = field->GetNumberOfComponent();
 639
 640             int nVectors = fop->GetNumberOfValues();
 641
 642             int block_size = ( nVectors / this->Internal->NumberOfPieces );
 643             int    start  = block_size * this->Internal->CurrentPieceNumber + 1;
 644             int    stride = block_size;
 645             int    count  = 1;
 646             int    blocksize = block_size;
 647             int    lastblocksize = (nVectors % this->Internal->NumberOfPieces);
 648             if ((this->Internal->NumberOfPieces ==
 649                  this->Internal->CurrentPieceNumber+1) && (lastblocksize != 0))
 650               {
 651               blocksize += lastblocksize;
 652               stride    += lastblocksize;
 653               }
 654             lastblocksize = 0;
 655
 656             vtkMedFilter *filter = vtkMedFilter::New();
 657             filter->SetFilterSizes( start, stride, count, blocksize, lastblocksize );
 658             fop->SetFilter(filter);
 659             }
 660           }
 661         }
 662       } // end IF
 663       }
 664     }
 665
 666 }
 667
 668 void  vtkMedReader::LinkMedInfo()
 669 {
 670   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 671       fieldfileit;
 672   // link the FieldOnProfile with the profiles
 673   fieldfileit = this->Internal->MedFiles.begin();
 674   while(fieldfileit != this->Internal->MedFiles.end())
 675     {
 676     vtkMedFile* fieldfile = fieldfileit->second;
 677     fieldfileit++;
 678
 679     for(int fid=0; fid<fieldfile->GetNumberOfField(); fid++)
 680       {
 681       vtkMedField* field = fieldfile->GetField(fid);
 682
 683       for(int sid = 0; sid< field->GetNumberOfFieldStep(); sid++)
 684         {
 685         vtkMedFieldStep* step = field->GetFieldStep(sid);
 686
 687         for(int foeid = 0; foeid < step->GetNumberOfFieldOverEntity(); foeid++)
 688           {
 689           vtkMedFieldOverEntity* fieldOverEntity = step->GetFieldOverEntity(foeid);
 690
 691           for(int fopid = 0; fopid < fieldOverEntity->GetNumberOfFieldOnProfile(); fopid++)
 692             {
 693             vtkMedFieldOnProfile* fop = fieldOverEntity->GetFieldOnProfile(fopid);
 694
 695             std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 696                 profilefileit = this->Internal->MedFiles.begin();
 697             while(profilefileit != this->Internal->MedFiles.end() && fop->GetProfile() == NULL)
 698               {
 699               vtkMedFile* profilefile = profilefileit->second;
 700               profilefileit++;
 701
 702               for(int pid = 0; pid < profilefile->GetNumberOfProfile(); pid++)
 703                 {
 704                 vtkMedProfile *profile = profilefile->GetProfile(pid);
 705                 if(strcmp(profile->GetName(), fop->GetProfileName()) == 0)
 706                   {
 707                   fop->SetProfile(profile);
 708                   break;
 709                   }
 710                 }
 711               }
 712             }
 713           }
 714         }
 715       }
 716     }
 717
 718   // first, add a familyOnEntityOnProfile to all FamilyOnEntity with a NULL
 719   // profile. This is used if no field is mapped to this support.
 720   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 721       meshfit = this->Internal->MedFiles.begin();
 722   while(meshfit != this->Internal->MedFiles.end())
 723     {
 724     vtkMedFile* meshfile = meshfit->second;
 725     meshfit++;
 726
 727     for(int mid=0; mid<meshfile->GetNumberOfMesh(); mid++)
 728       {
 729       vtkMedMesh* mesh = meshfile->GetMesh(mid);
 730
 731       for(int gid=0; gid<mesh->GetNumberOfGridStep(); gid++)
 732         {
 733         vtkMedGrid* grid = mesh->GetGridStep(gid);
 734         // read point family data and create EntityArrays
 735
 736         for(int eid=0; eid < grid->GetNumberOfEntityArray(); eid++)
 737           {
 738           vtkMedEntityArray* array = grid->GetEntityArray(eid);
 739
 740           for(int fid=0; fid < array->GetNumberOfFamilyOnEntity(); fid++)
 741             {
 742             vtkMedFamilyOnEntity* foe = array->GetFamilyOnEntity(fid);
 743             if(foe->GetFamilyOnEntityOnProfile((vtkMedProfile*)NULL) == NULL)
 744               {
 745               vtkMedFamilyOnEntityOnProfile* foep =
 746                   vtkMedFamilyOnEntityOnProfile::New();
 747               foep->SetFamilyOnEntity(foe);
 748               foep->SetProfile(NULL);
 749               foe->AddFamilyOnEntityOnProfile(foep);
 750               foep->Delete();
 751               }
 752             }//family on entity
 753           }//entity array
 754         }// grid step
 755       }//mesh
 756     }//mesh file
 757
 758   fieldfileit = this->Internal->MedFiles.begin();
 759   while(fieldfileit != this->Internal->MedFiles.end())
 760     {
 761     vtkMedFile* fieldfile = fieldfileit->second;
 762     fieldfileit++;
 763
 764     for(int fieldid=0; fieldid < fieldfile->GetNumberOfField(); fieldid++)
 765       {
 766       vtkMedField* field = fieldfile->GetField(fieldid);
 767
 768       for(int fstepid=0; fstepid < field->GetNumberOfFieldStep(); fstepid++)
 769         {
 770         vtkMedFieldStep* step = field->GetFieldStep(fstepid);
 771
 772         vtkMedComputeStep meshcs = step->GetMeshComputeStep();
 773
 774         for(int foeid=0; foeid<step->GetNumberOfFieldOverEntity() ;foeid++)
 775           {
 776           vtkMedFieldOverEntity* fieldOverEntity = step->GetFieldOverEntity(foeid);
 777           const vtkMedEntity& fieldentity = fieldOverEntity->GetEntity();
 778
 779           for (int fopid = 0;
 780                fopid < fieldOverEntity->GetNumberOfFieldOnProfile(); fopid++)
 781             {
 782             vtkMedFieldOnProfile* fop =
 783                 fieldOverEntity->GetFieldOnProfile(fopid);
 784
 785             std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 786                 meshfileit = this->Internal->MedFiles.begin();
 787             while(meshfileit != this->Internal->MedFiles.end())
 788               {
 789               vtkMedFile* meshfile = meshfileit->second;
 790               meshfileit++;
 791
 792               if(field->GetLocal() == 1 && (meshfile != fieldfile))
 793                 continue;
 794
 795               for(int mid=0; mid<meshfile->GetNumberOfMesh(); mid++)
 796                 {
 797                 vtkMedMesh* mesh = meshfile->GetMesh(mid);
 798
 799                 // the field must be on this mesh.
 800                 if(strcmp(mesh->GetName(),
 801                           field->GetMeshName()) != 0)
 802                   continue;
 803
 804                 vtkMedGrid* grid = mesh->GetGridStep(meshcs);
 805                 if(grid == NULL)
 806                   {
 807                   vtkErrorMacro("the field " << field->GetName()
 808                                 << " at step iteration:"
 809                                 << step->GetComputeStep().IterationIt
 810                                 << " and time "
 811                                 << step->GetComputeStep().TimeIt
 812                                 << " uses mesh at step "
 813                                 << meshcs.TimeIt << " " << meshcs.IterationIt
 814                                 << "which does not exists!");
 815                   continue;
 816                   }
 817
 818                 for(int eid=0; eid < grid->GetNumberOfEntityArray(); eid++)
 819                   {
 820                   vtkMedEntityArray* array = grid->GetEntityArray(eid);
 821
 822                   // if the support is on points,
 823                   // the field must also be on points
 824                   if(array->GetEntity().EntityType == MED_NODE &&
 825                      fieldentity.EntityType != MED_NODE)
 826                     continue;
 827
 828                   if(array->GetEntity().EntityType != MED_NODE &&
 829                      fieldentity.EntityType == MED_NODE)
 830                     continue;
 831
 832                   // for fields not on points, the geometry type
 833                   // of the support must match
 834                   if(array->GetEntity().EntityType != MED_NODE &&
 835                      array->GetEntity().GeometryType != fieldentity.GeometryType)
 836                     continue;
 837
 838                   for(int fid = 0; fid < array->GetNumberOfFamilyOnEntity(); fid++)
 839                     {
 840                     vtkMedFamilyOnEntity* foe = array->GetFamilyOnEntity(fid);
 841                     if(foe->GetFamilyOnEntityOnProfile(fop->GetProfile()) == NULL)
 842                       {
 843                       vtkMedFamilyOnEntityOnProfile* foep =
 844                           vtkMedFamilyOnEntityOnProfile::New();
 845                       foep->SetProfile(fop->GetProfile());
 846                       foep->SetFamilyOnEntity(foe);
 847                       foe->AddFamilyOnEntityOnProfile(foep);
 848                       foep->Delete();
 849                       }
 850                     // also add the family on entity with no profile.
 851                     if(foe->GetFamilyOnEntityOnProfile((vtkMedProfile*)NULL) == NULL)
 852                       {
 853                       vtkMedFamilyOnEntityOnProfile* foep =
 854                           vtkMedFamilyOnEntityOnProfile::New();
 855                       foep->SetProfile(NULL);
 856                       foep->SetFamilyOnEntity(foe);
 857                       foe->AddFamilyOnEntityOnProfile(foep);
 858                       foep->Delete();
 859                       }
 860                     }//familyOnEntity
 861                   }//entityArray
 862                 }//mesh
 863               }//mesh file
 864             }//field on profile
 865           }//fieldOverEntity
 866         }//fiedstep
 867       }// fieldid
 868     }//fieldfileit
 869
 870   // Now, link localizations and interpolations
 871   fieldfileit = this->Internal->MedFiles.begin();
 872   while(fieldfileit != this->Internal->MedFiles.end())
 873     {
 874     vtkMedFile* fieldfile = fieldfileit->second;
 875     fieldfileit++;
 876
 877     for(int locid = 0; locid < fieldfile->GetNumberOfLocalization(); locid++)
 878       {
 879       vtkMedLocalization* loc = fieldfile->GetLocalization(locid);
 880
 881       for(int fid = 0; fid < fieldfile->GetNumberOfField() &&
 882                     loc->GetInterpolation() == NULL; fid++)
 883         {
 884         vtkMedField* field = fieldfile->GetField(fid);
 885         for(int interpid = 0; interpid < field->GetNumberOfInterpolation();
 886         interpid++)
 887           {
 888           vtkMedInterpolation* interp = field->GetInterpolation(interpid);
 889           if(strcmp(loc->GetInterpolationName(),
 890                     interp->GetName()) == 0)
 891             {
 892             loc->SetInterpolation(interp);
 893             break;
 894             }
 895           }
 896         }
 897       }
 898     }
 899
 900   // now that the interpolation is set, build the  shape functions.
 901   fieldfileit = this->Internal->MedFiles.begin();
 902   while(fieldfileit != this->Internal->MedFiles.end())
 903     {
 904     vtkMedFile* fieldfile = fieldfileit->second;
 905     fieldfileit++;
 906
 907     for(int locid = 0; locid < fieldfile->GetNumberOfLocalization(); locid++)
 908       {
 909       vtkMedLocalization* loc = fieldfile->GetLocalization(locid);
 910       loc->BuildShapeFunction();
 911       }
 912     }
 913
 914   // set the supportmesh pointer in the structural element
 915   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
 916       fileit = this->Internal->MedFiles.begin();
 917   while(fileit != this->Internal->MedFiles.end())
 918     {
 919     vtkMedFile* file = fileit->second;
 920     fileit++;
 921
 922     for(int structelemit = 0;
 923         structelemit < file->GetNumberOfStructElement();
 924         structelemit++)
 925       {
 926       vtkMedStructElement* structElem =
 927           file->GetStructElement(structelemit);
 928
 929       for(int supportmeshit = 0;
 930           supportmeshit < file->GetNumberOfSupportMesh();
 931           supportmeshit++)
 932         {
 933         vtkMedMesh* supportMesh =
 934             file->GetSupportMesh(supportmeshit);
 935
 936         if(strcmp(supportMesh->GetName(), structElem->GetName()) == 0 )
 937           {
 938           structElem->SetSupportMesh(supportMesh);
 939           break;
 940           }
 941         }
 942       }
 943     }
 944
 945   // set the pointer to the profile used by the constant attributes
 946   fileit = this->Internal->MedFiles.begin();
 947   while(fileit != this->Internal->MedFiles.end())
 948   {
 949   vtkMedFile* file = fileit->second;
 950   fileit++;
 951
 952   for(int structelemit = 0;
 953       structelemit < file->GetNumberOfStructElement();
 954       structelemit++)
 955     {
 956     vtkMedStructElement* structElem =
 957         file->GetStructElement(structelemit);
 958
 959     for(int cstattit = 0; cstattit < structElem->GetNumberOfConstantAttribute(); cstattit++)
 960       {
 961       vtkMedConstantAttribute* cstatt = structElem->GetConstantAttribute(cstattit);
 962
 963       for(int profit = 0;
 964           profit < file->GetNumberOfProfile();
 965           profit++)
 966         {
 967         vtkMedProfile* profile =
 968             file->GetProfile(profit);
 969
 970         if(strcmp(profile->GetName(), cstatt->GetProfileName()) == 0 )
 971           {
 972           cstatt->SetProfile(profile);
 973           break;
 974           }
 975         }
 976       }
 977     }
 978   }
 979
 980   meshfit = this->Internal->MedFiles.begin();
 981   while(meshfit != this->Internal->MedFiles.end())
 982   {
 983   vtkMedFile* meshfile = meshfit->second;
 984   meshfit++;
 985
 986   for(int mid=0; mid<meshfile->GetNumberOfMesh(); mid++)
 987     {
 988     vtkMedMesh* mesh = meshfile->GetMesh(mid);
 989
 990     for(int gid=0; gid<mesh->GetNumberOfGridStep(); gid++)
 991       {
 992       vtkMedGrid* grid = mesh->GetGridStep(gid);
 993       // read point family data and create EntityArrays
 994
 995       for(int eid=0; eid < grid->GetNumberOfEntityArray(); eid++)
 996         {
 997         vtkMedEntityArray* array = grid->GetEntityArray(eid);
 998         if(array->GetEntity().EntityType != MED_STRUCT_ELEMENT)
 999           continue;
1000
1001         for(int structelemit = 0; structelemit < meshfile->GetNumberOfStructElement(); structelemit++)
1002           {
1003           vtkMedStructElement* structelem = meshfile->GetStructElement(structelemit);
1004           if(structelem->GetGeometryType() == array->GetEntity().GeometryType)
1005             {
1006             array->SetStructElement(structelem);
1007             break;
1008             }
1009           }
1010         }
1011       }
1012     }
1013   }
1014 }
1015
1016 int vtkMedReader::GetFrequencyArrayStatus(const char* name)
1017 {
1018   return this->Frequencies->GetKeyStatus(name);
1019 }
1020
1021 void vtkMedReader::SetFrequencyArrayStatus(const char* name, int status)
1022 {
1023   if(this->Frequencies->GetKeyStatus(name) == status)
1024     {
1025     return;
1026     }
1027
1028   this->Frequencies->SetKeyStatus(name, status);
1029
1030   this->Modified();
1031 }
1032
1033 int vtkMedReader::GetNumberOfFrequencyArrays()
1034 {
1035   return this->Frequencies->GetNumberOfKey();
1036 }
1037
1038 const char* vtkMedReader::GetFrequencyArrayName(int index)
1039 {
1040   return this->Frequencies->GetKey(index);
1041 }
1042
1043 struct compTimes
1044 {
1045   bool operator()(pair<double, int> i, pair<double, int> j)
1046   {
1047     if(i.first!=j.first)
1048       return (i.first<j.first);
1049     return i.second<j.second;
1050   }
1051 };
1052
1053 vtkDoubleArray* vtkMedReader::GetAvailableTimes()
1054 {
1055   this->AvailableTimes->Initialize();
1056   this->AvailableTimes->SetNumberOfComponents(1);
1057
1058   std::set<std::string> newFrequencies;
1059
1060   int tid = 0;
1061   std::map<med_float, std::set<med_int> >::iterator it =
1062       this->Internal->GlobalComputeStep.begin();
1063   while(it != this->Internal->GlobalComputeStep.end())
1064     {
1065     double time = it->first;
1066     this->AvailableTimes->InsertNextValue(time);
1067     string name = vtkMedUtilities::GetModeKey(tid, time, this->Internal->GlobalComputeStep.size()-1);
1068     this->Frequencies->AddKey(name.c_str());
1069     newFrequencies.insert(name);
1070     tid++;
1071     it++;
1072     }
1073
1074   // now check if old frequencies have been removed
1075   for(int f = 0; f < this->Frequencies->GetNumberOfKey(); f++)
1076     {
1077     const char* name = this->Frequencies->GetKey(f);
1078     if(newFrequencies.find(name) == newFrequencies.end())
1079       {
1080       this->Frequencies->RemoveKeyByIndex(f);
1081       f--;
1082       }
1083     }
1084
1085   return this->AvailableTimes;
1086 }
1087
1088 void vtkMedReader::ChooseRealAnimationMode()
1089 {
1090   if(this->AnimationMode!=Default)
1091     {
1092     this->Internal->RealAnimationMode=this->AnimationMode;
1093     return;
1094     }
1095
1096   // if there is exactly one physical time and more than one iteration
1097   // set the animation mode to iteration, else default to physical time.
1098   if (this->Internal->GlobalComputeStep.size() == 1 &&
1099       this->Internal->GlobalComputeStep[0].size() > 1)
1100     {
1101     this->Internal->RealAnimationMode=Iteration;
1102     return;
1103     }
1104
1105   this->Internal->RealAnimationMode=PhysicalTime;
1106 }
1107
1108 int vtkMedReader::GetEntityStatus(const vtkMedEntity& entity)
1109 {
1110   if (entity.EntityType==MED_NODE)
1111     return 1;
1112   if(entity.EntityType == MED_DESCENDING_FACE
1113      || entity.EntityType == MED_DESCENDING_EDGE)
1114     return 0;
1115
1116   return this->Entities->GetKeyStatus(vtkMedUtilities::EntityKey(entity).c_str());
1117 }
1118
1119 int vtkMedReader::GetFamilyStatus(vtkMedMesh* mesh, vtkMedFamily* family)
1120 {
1121   if(!mesh||!family)
1122     return 0;
1123
1124   if(this->Internal->GroupSelectionMTime > this->Internal->FamilySelectionMTime)
1125     {
1126     this->SelectFamiliesFromGroups();
1127     }
1128
1129   int status =  this->Internal->Families->GetKeyStatus(vtkMedUtilities::FamilyKey(
1130       mesh->GetName(), family->GetPointOrCell(),
1131       family->GetName()).c_str());
1132
1133   return status;
1134 }
1135
1136 int vtkMedReader::IsMeshSelected(vtkMedMesh* mesh)
1137 {
1138   for(int fam=0; fam<mesh->GetNumberOfPointFamily(); fam++)
1139     {
1140     if(this->GetFamilyStatus(mesh, mesh->GetPointFamily(fam))!=0)
1141       return 1;
1142     }
1143
1144   for(int fam=0; fam<mesh->GetNumberOfCellFamily(); fam++)
1145     {
1146     if(this->GetFamilyStatus(mesh, mesh->GetCellFamily(fam))!=0)
1147       return 1;
1148     }
1149   return 0;
1150 }
1151
1152 void vtkMedReader::GatherComputeSteps()
1153 {
1154   this->Internal->GlobalComputeStep.clear();
1155
1156   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1157       fieldfileit = this->Internal->MedFiles.begin();
1158   while(fieldfileit != this->Internal->MedFiles.end())
1159     {
1160     vtkMedFile* file = fieldfileit->second;
1161     fieldfileit++;
1162
1163     // first loop over all fields to gather their compute steps
1164     for(int fieldId=0; fieldId<file->GetNumberOfField(); fieldId++)
1165       {
1166       vtkMedField* field=file->GetField(fieldId);
1167
1168       for(int stepId=0; stepId<field->GetNumberOfFieldStep(); stepId++)
1169         {
1170         vtkMedFieldStep* step=field->GetFieldStep(stepId);
1171         const vtkMedComputeStep& cs = step->GetComputeStep();
1172         this->Internal->GlobalComputeStep[cs.TimeOrFrequency].insert(cs.IterationIt);
1173         }
1174       }//fields
1175
1176     // then loop over all meshes to gather their grid steps too.
1177     // for meshes, do not add the MED_UNDEF_DT time
1178     for(int meshId=0; meshId<file->GetNumberOfMesh(); meshId++)
1179       {
1180       vtkMedMesh* mesh=file->GetMesh(meshId);
1181
1182       for(int stepId=0; stepId<mesh->GetNumberOfGridStep(); stepId++)
1183         {
1184         vtkMedGrid* grid=mesh->GetGridStep(stepId);
1185         const vtkMedComputeStep& cs = grid->GetComputeStep();
1186         if(cs.TimeOrFrequency != MED_UNDEF_DT || cs.TimeIt != MED_NO_DT)
1187           {
1188           this->Internal->GlobalComputeStep[cs.TimeOrFrequency].insert(cs.IterationIt);
1189           }
1190         }
1191       }//mesh
1192     }
1193   if(this->Internal->GlobalComputeStep.size() == 0)
1194     {
1195     this->Internal->GlobalComputeStep[MED_UNDEF_DT].insert(MED_NO_IT);
1196     }
1197 }
1198
1199 int vtkMedReader::IsFieldSelected(vtkMedField* field)
1200 {
1201   return this->IsPointFieldSelected(field)||this->IsCellFieldSelected(field)
1202       ||this->IsQuadratureFieldSelected(field) || this->IsElnoFieldSelected(field);
1203 }
1204
1205 int vtkMedReader::IsPointFieldSelected(vtkMedField* field)
1206 {
1207   return field->GetFieldType()==vtkMedField::PointField
1208       &&this->GetPointFieldArrayStatus(vtkMedUtilities::SimplifyName(
1209           field->GetName()).c_str());
1210 }
1211
1212 int vtkMedReader::IsCellFieldSelected(vtkMedField* field)
1213 {
1214   return field->GetFieldType()==vtkMedField::CellField
1215       &&this->GetCellFieldArrayStatus(vtkMedUtilities::SimplifyName(
1216           field->GetName()).c_str());
1217 }
1218
1219 vtkMedProfile* vtkMedReader::GetProfile(const char* pname)
1220 {
1221   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1222       fileit = this->Internal->MedFiles.begin();
1223   while(fileit != this->Internal->MedFiles.end())
1224     {
1225     vtkMedFile* file = fileit->second;
1226     fileit++;
1227     vtkMedProfile* profile = file->GetProfile(pname);
1228     if(profile != NULL)
1229       return profile;
1230     }
1231   return NULL;
1232 }
1233
1234 vtkMedLocalization* vtkMedReader::GetLocalization(const char* lname)
1235 {
1236   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1237       fileit = this->Internal->MedFiles.begin();
1238   while(fileit != this->Internal->MedFiles.end())
1239     {
1240     vtkMedFile* file = fileit->second;
1241     fileit++;
1242     vtkMedLocalization* loc = file->GetLocalization(lname);
1243     if(loc != NULL)
1244       return loc;
1245     }
1246   return NULL;
1247
1248 }
1249
1250 int vtkMedReader::GetLocalizationKey(vtkMedFieldOnProfile* fop)
1251 {
1252   vtkMedLocalization* def=this->GetLocalization(fop->GetLocalizationName());
1253
1254   // This is not a quadrature field with explicit definition.
1255   // There are two possible cases : either the intergration point is
1256   // at the center of the cell
1257   //1 quadrature point at the cell center
1258   int nloc = 0;
1259   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1260       fileit = this->Internal->MedFiles.begin();
1261   while(fileit != this->Internal->MedFiles.end())
1262     {
1263     vtkMedFile* file = fileit->second;
1264     fileit++;
1265
1266     if(def && def->GetParentFile() == file)
1267       return nloc + def->GetMedIterator() - 1;
1268
1269     nloc += file->GetNumberOfLocalization();
1270     }
1271
1272   // center of a cell
1273   if(fop->GetNumberOfIntegrationPoint()==1)
1274     return nloc + 1 + fop->GetParentFieldOverEntity()->GetEntity().GeometryType;
1275
1276   // or it is an elno field (field stored on nodes of the cells,
1277   // but with discontinuities at the vertices)
1278   return -fop->GetParentFieldOverEntity()->GetEntity().GeometryType;//ELNO
1279 }
1280
1281 int vtkMedReader::IsQuadratureFieldSelected(vtkMedField* field)
1282 {
1283   return field->GetFieldType()==vtkMedField::QuadratureField
1284       &&this->GetQuadratureFieldArrayStatus(vtkMedUtilities::SimplifyName(
1285           field->GetName()).c_str());
1286 }
1287
1288 int vtkMedReader::IsElnoFieldSelected(vtkMedField* field)
1289 {
1290   return field->GetFieldType()==vtkMedField::ElnoField
1291       &&this->GetElnoFieldArrayStatus(vtkMedUtilities::SimplifyName(
1292           field->GetName()).c_str());
1293 }
1294
1295 // Description:
1296 // Give the animation steps to the pipeline
1297 void vtkMedReader::AdvertiseTime(vtkInformation* info)
1298 {
1299   this->ChooseRealAnimationMode();
1300
1301   if(this->Internal->RealAnimationMode==PhysicalTime)
1302     {
1303     // I advertise the union of all times available
1304     // in all selected fields and meshes
1305     set<double> timeset;
1306
1307     std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1308         fieldfileit = this->Internal->MedFiles.begin();
1309     while(fieldfileit != this->Internal->MedFiles.end())
1310       {
1311       vtkMedFile* file = fieldfileit->second;
1312       fieldfileit++;
1313
1314       // first loop over all fields to gather their compute steps
1315       for(int fieldId=0; fieldId<file->GetNumberOfField(); fieldId++)
1316         {
1317         vtkMedField* field=file->GetField(fieldId);
1318
1319         if(!this->IsFieldSelected(field))
1320           continue;
1321
1322         field->GatherFieldTimes(timeset);
1323         }//fields
1324
1325       // then loop over all meshes to gather their grid steps too.
1326       for(int meshId=0; meshId<file->GetNumberOfMesh(); meshId++)
1327         {
1328         vtkMedMesh* mesh=file->GetMesh(meshId);
1329
1330         if(!this->IsMeshSelected(mesh))
1331           continue;
1332
1333         mesh->GatherGridTimes(timeset);
1334         }//meshes
1335       }
1336
1337     if(timeset.size() > 0)
1338       {
1339       // remove MED_UNDEF_DT if there are other time step
1340       if(timeset.size() > 1)
1341         timeset.erase(MED_UNDEF_DT);
1342
1343       vector<double> times;
1344       set<double>::iterator it = timeset.begin();
1345       while(it != timeset.end())
1346         {
1347         times.push_back(*it);
1348         it++;
1349         }
1350       sort(times.begin(), times.end());
1351
1352       info->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &times[0],
1353           times.size());
1354       double timeRange[2];
1355       timeRange[0]=times[0];
1356       timeRange[1]=times[times.size()-1];
1357       info->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), &timeRange[0],
1358           2);
1359       }
1360     else
1361       {
1362       info->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
1363       info->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
1364       }
1365     }
1366   else if(this->Internal->RealAnimationMode==Iteration)
1367     {
1368     // I advertise the union of all iterations available at the given
1369     // Time for all selected fields.
1370     set<med_int> iterationsets;
1371     med_float time = MED_UNDEF_DT;
1372     if(this->TimeIndexForIterations >= 0 &&
1373        this->TimeIndexForIterations <
1374        this->AvailableTimes->GetNumberOfTuples())
1375       {
1376       time = this->AvailableTimes->
1377                      GetValue((vtkIdType)this->TimeIndexForIterations);
1378       }
1379
1380     std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1381         fieldfileit = this->Internal->MedFiles.begin();
1382     while(fieldfileit != this->Internal->MedFiles.end())
1383       {
1384       vtkMedFile* file = fieldfileit->second;
1385       fieldfileit++;
1386
1387       for(int fieldId=0; fieldId<file->GetNumberOfField(); fieldId++)
1388         {
1389         vtkMedField* field=file->GetField(fieldId);
1390         if(!this->IsFieldSelected(field))
1391           continue;
1392
1393         field->GatherFieldIterations(time, iterationsets);
1394         }
1395       // then loop over all meshes to gather their grid steps too.
1396       for(int meshId=0; meshId<file->GetNumberOfMesh(); meshId++)
1397         {
1398         vtkMedMesh* mesh=file->GetMesh(meshId);
1399
1400         if(!this->IsMeshSelected(mesh))
1401           continue;
1402
1403         mesh->GatherGridIterations(time, iterationsets);
1404         }//meshes
1405       }
1406
1407     if(iterationsets.size()>0)
1408       {
1409       // remove MED_NO_IT if there are other available iterations.
1410       if(iterationsets.size()>1)
1411         iterationsets.erase(MED_NO_IT);
1412
1413       vector<double> iterations;
1414       set<med_int>::iterator it=iterationsets.begin();
1415       while(it!=iterationsets.end())
1416         {
1417         iterations.push_back((double)*it);
1418         it++;
1419         }
1420       sort(iterations.begin(), iterations.end());
1421       info->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &iterations[0],
1422           iterations.size());
1423       double timeRange[2];
1424       timeRange[0]=iterations[0];
1425       timeRange[1]=iterations[iterations.size()-1];
1426       info->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), &timeRange[0],
1427           2);
1428       }
1429     else
1430       {
1431       info->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
1432       info->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
1433       }
1434     }
1435   else
1436     {
1437     info->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
1438     info->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
1439     }
1440 }
1441
1442 vtkIdType vtkMedReader::GetFrequencyIndex(double freq)
1443 {
1444   return this->AvailableTimes->LookupValue(freq);
1445 }
1446
1447 int vtkMedReader::RequestDataObject(vtkInformation* request,
1448     vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
1449 {
1450   vtkInformation *info = outputVector->GetInformationObject(0);
1451   if (vtkMultiBlockDataSet::SafeDownCast(
1452       info->Get(vtkDataObject::DATA_OBJECT())))
1453     {
1454     // The output is already created
1455     return 1;
1456     }
1457   else
1458     {
1459     vtkMultiBlockDataSet* output=vtkMultiBlockDataSet::New();
1460     this->GetExecutive()->SetOutputData(0, output);
1461     output->Delete();
1462     this->GetOutputPortInformation(0)->Set(vtkDataObject::DATA_EXTENT_TYPE(),
1463         output->GetExtentType());
1464     }
1465   return 1;
1466 }
1467
1468 void vtkMedReader::ClearSelections()
1469 {
1470   this->PointFields->Initialize();
1471   this->CellFields->Initialize();
1472   this->QuadratureFields->Initialize();
1473   this->ElnoFields->Initialize();
1474
1475   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1476       fieldfileit = this->Internal->MedFiles.begin();
1477   while(fieldfileit != this->Internal->MedFiles.end())
1478     {
1479     vtkMedFile* file = fieldfileit->second;
1480     fieldfileit++;
1481
1482     for(int index=0; index < file->GetNumberOfField(); index++)
1483       {
1484       vtkMedField* field = file->GetField(index);
1485       switch(field->GetFieldType())
1486         {
1487         case vtkMedField::PointField :
1488         this->PointFields->AddKey(vtkMedUtilities::SimplifyName(
1489               field->GetName()).c_str());
1490         break;
1491         case vtkMedField::CellField :
1492         this->CellFields->AddKey(vtkMedUtilities::SimplifyName(
1493               field->GetName()).c_str());
1494         break;
1495         case vtkMedField::QuadratureField :
1496         this->QuadratureFields->AddKey(vtkMedUtilities::SimplifyName(
1497               field->GetName()).c_str());
1498         break;
1499         case vtkMedField::ElnoField :
1500         this->ElnoFields->AddKey(vtkMedUtilities::SimplifyName(
1501               field->GetName()).c_str());
1502         break;
1503         }
1504       }
1505
1506     this->Internal->Families->Initialize();
1507     this->Groups->Initialize();
1508     for(int meshIndex=0; meshIndex < file->GetNumberOfMesh(); meshIndex++)
1509       {
1510       vtkMedMesh* mesh = file->GetMesh(meshIndex);
1511       for(int famIndex=0; famIndex<mesh->GetNumberOfPointFamily(); famIndex++)
1512         {
1513         vtkMedFamily* fam=mesh->GetPointFamily(famIndex);
1514
1515         int ng=fam->GetNumberOfGroup();
1516         for(int gindex=0; gindex<ng; gindex++)
1517           {
1518           vtkMedGroup* group=fam->GetGroup(gindex);
1519           string gname=vtkMedUtilities::GroupKey(mesh->GetName(),
1520               fam->GetPointOrCell(), group->GetName());
1521
1522           this->Groups->AddKey(gname.c_str());
1523           this->Groups->SetKeyStatus(gname.c_str(), 0);
1524           }
1525         }
1526       for(int famIndex=0; famIndex<mesh->GetNumberOfCellFamily(); famIndex++)
1527         {
1528         vtkMedFamily* fam=mesh->GetCellFamily(famIndex);
1529
1530         int ng=fam->GetNumberOfGroup();
1531         for(int gindex=0; gindex<ng; gindex++)
1532           {
1533           vtkMedGroup* group=fam->GetGroup(gindex);
1534           string gname=vtkMedUtilities::GroupKey(mesh->GetName(),
1535               fam->GetPointOrCell(), group->GetName());
1536
1537           this->Groups->AddKey(gname.c_str());
1538           this->Groups->SetKeyStatus(gname.c_str(), 1);
1539           }
1540         }
1541       }
1542     this->Internal->GroupSelectionMTime.Modified();
1543
1544     for(int meshIndex=0; meshIndex< file->GetNumberOfMesh(); meshIndex++)
1545       {
1546       if(file->GetMesh(meshIndex)->GetNumberOfGridStep() == 0)
1547         continue;
1548
1549       vtkMedGrid* grid=file->GetMesh(meshIndex)->GetGridStep(0);
1550
1551       for(int entityIndex=0; entityIndex<grid->GetNumberOfEntityArray();
1552         entityIndex++)
1553         {
1554         vtkMedEntityArray* array=grid->GetEntityArray(entityIndex);
1555         string name=vtkMedUtilities::EntityKey(array->GetEntity());
1556         this->Entities->AddKey(name.c_str());
1557         }
1558       }
1559     }
1560   this->Modified();
1561 }
1562
1563 void vtkMedReader::SelectFamiliesFromGroups()
1564 {
1565   this->Internal->Families->Initialize();
1566
1567   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1568       meshfileit = this->Internal->MedFiles.begin();
1569   while(meshfileit != this->Internal->MedFiles.end())
1570     {
1571     vtkMedFile* file = meshfileit->second;
1572     meshfileit++;
1573
1574     for(int meshIndex=0; meshIndex < file->GetNumberOfMesh(); meshIndex++)
1575       {
1576       vtkMedMesh* mesh = file->GetMesh(meshIndex);
1577       for(int famIndex=0; famIndex<mesh->GetNumberOfFamily(); famIndex++)
1578         {
1579         vtkMedFamily* fam=mesh->GetFamily(famIndex);
1580         string name=vtkMedUtilities::FamilyKey(mesh->GetName(),
1581             fam->GetPointOrCell(), fam->GetName());
1582
1583         this->Internal->Families->SetKeyStatus(name.c_str(), 0);
1584
1585         for(int gindex=0; gindex<fam->GetNumberOfGroup(); gindex++)
1586           {
1587           vtkMedGroup* group=fam->GetGroup(gindex);
1588           string gname=vtkMedUtilities::GroupKey(mesh->GetName(),
1589               fam->GetPointOrCell(), group->GetName());
1590           int state=this->Groups->GetKeyStatus(gname.c_str());
1591
1592           if(state)
1593             {
1594             this->SetFamilyStatus(name.c_str(), 1);
1595             }
1596           }
1597         }
1598       }
1599   }
1600
1601   this->Internal->FamilySelectionMTime.Modified();
1602 }
1603
1604 int vtkMedReader::GetNumberOfPointFieldArrays()
1605 {
1606   return this->PointFields->GetNumberOfKey();
1607 }
1608
1609 const char*
1610 vtkMedReader::GetPointFieldArrayName(int index)
1611 {
1612   return this->PointFields->GetKey(index);
1613 }
1614
1615 int vtkMedReader::GetPointFieldArrayStatus(const char* name)
1616 {
1617   return this->PointFields->GetKeyStatus(name);
1618 }
1619
1620 void vtkMedReader::SetPointFieldArrayStatus(const char* name, int status)
1621 {
1622   if(this->PointFields->KeyExists(name)&&this->PointFields->GetKeyStatus(
1623       name)==status)
1624     return;
1625
1626   this->PointFields->SetKeyStatus(name, status);
1627
1628   this->Modified();
1629 }
1630
1631 int vtkMedReader::GetNumberOfCellFieldArrays()
1632 {
1633   return this->CellFields->GetNumberOfKey();
1634 }
1635
1636 const char*
1637 vtkMedReader::GetCellFieldArrayName(int index)
1638 {
1639   return this->CellFields->GetKey(index);
1640 }
1641
1642 int vtkMedReader::GetCellFieldArrayStatus(const char* name)
1643 {
1644   return this->CellFields->GetKeyStatus(name);
1645 }
1646
1647 void vtkMedReader::SetCellFieldArrayStatus(const char* name, int status)
1648 {
1649   if(this->CellFields->KeyExists(name)&&this->CellFields->GetKeyStatus(
1650       name)==status)
1651     return;
1652
1653   this->CellFields->SetKeyStatus(name, status);
1654
1655   this->Modified();
1656 }
1657
1658 int vtkMedReader::GetNumberOfQuadratureFieldArrays()
1659 {
1660   return this->QuadratureFields->GetNumberOfKey();
1661 }
1662
1663 const char* vtkMedReader::GetQuadratureFieldArrayName(int index)
1664 {
1665   return this->QuadratureFields->GetKey(index);
1666 }
1667
1668 int vtkMedReader::GetQuadratureFieldArrayStatus(const char* name)
1669 {
1670   return this->QuadratureFields->GetKeyStatus(name);
1671 }
1672
1673 void vtkMedReader::SetQuadratureFieldArrayStatus(const char* name, int status)
1674 {
1675   if(this->QuadratureFields->KeyExists(name)
1676       &&this->QuadratureFields->GetKeyStatus(name)==status)
1677     return;
1678
1679   this->QuadratureFields->SetKeyStatus(name, status);
1680
1681   this->Modified();
1682 }
1683
1684 int vtkMedReader::GetNumberOfElnoFieldArrays()
1685 {
1686   return this->ElnoFields->GetNumberOfKey();
1687 }
1688
1689 const char* vtkMedReader::GetElnoFieldArrayName(int index)
1690 {
1691   return this->ElnoFields->GetKey(index);
1692 }
1693
1694 int vtkMedReader::GetElnoFieldArrayStatus(const char* name)
1695 {
1696   return this->ElnoFields->GetKeyStatus(name);
1697 }
1698
1699 void vtkMedReader::SetElnoFieldArrayStatus(const char* name, int status)
1700 {
1701   if(this->ElnoFields->KeyExists(name)
1702       &&this->ElnoFields->GetKeyStatus(name)==status)
1703     return;
1704
1705   this->ElnoFields->SetKeyStatus(name, status);
1706
1707   this->Modified();
1708 }
1709
1710 void vtkMedReader::SetEntityStatus(const char* name, int status)
1711 {
1712   if(this->Entities->KeyExists(name)&&this->Entities->GetKeyStatus(name)
1713       ==status)
1714     return;
1715
1716   this->Entities->SetKeyStatus(name, status);
1717
1718   this->Modified();
1719 }
1720
1721 void vtkMedReader::SetFamilyStatus(const char* name, int status)
1722 {
1723   if(this->Internal->Families->KeyExists(name)
1724       &&this->Internal->Families->GetKeyStatus(name)==status)
1725     return;
1726
1727   this->Internal->Families->SetKeyStatus(name, status);
1728 }
1729
1730 void vtkMedReader::SetGroupStatus(const char* name, int status)
1731 {
1732
1733   if(this->Groups->KeyExists(name)&&this->Groups->GetKeyStatus(name)
1734       ==status)
1735     return;
1736
1737   this->Groups->SetKeyStatus(name, status);
1738
1739   this->Modified();
1740
1741   this->Internal->GroupSelectionMTime.Modified();
1742 }
1743
1744 int vtkMedReader::GetGroupStatus(const char* key)
1745 {
1746   return this->Groups->GetKeyStatus(key);
1747 }
1748
1749 void vtkMedReader::AddQuadratureSchemeDefinition(vtkInformation* info,
1750     vtkMedLocalization* loc)
1751 {
1752   if(info==NULL||loc==NULL)
1753     return;
1754
1755   vtkInformationQuadratureSchemeDefinitionVectorKey *key=
1756       vtkQuadratureSchemeDefinition::DICTIONARY();
1757
1758   vtkQuadratureSchemeDefinition* def=vtkQuadratureSchemeDefinition::New();
1759   int cellType=vtkMedUtilities::GetVTKCellType(loc->GetGeometryType());
1760   def->Initialize(cellType, vtkMedUtilities::GetNumberOfPoint(
1761       loc->GetGeometryType()), loc->GetNumberOfQuadraturePoint(),
1762       (double*)loc->GetShapeFunction()->GetVoidPointer(0),
1763       (double*)loc->GetWeights()->GetVoidPointer(0));
1764   key->Set(info, def, cellType);
1765   def->Delete();
1766
1767 }
1768
1769 void vtkMedReader::LoadConnectivity(vtkMedEntityArray* array)
1770 {
1771   vtkMedGrid* grid = array->GetParentGrid();
1772   array->LoadConnectivity();
1773   if (array->GetConnectivity()==MED_NODAL||vtkMedUtilities::GetDimension(
1774       array->GetEntity().GeometryType)<2
1775       || grid->GetParentMesh()->GetMeshType() == MED_STRUCTURED_MESH)
1776     return;
1777
1778   vtkMedEntity subentity;
1779   subentity.EntityType = vtkMedUtilities::GetSubType(array->GetEntity().EntityType);
1780
1781   vtkMedUnstructuredGrid* ugrid = vtkMedUnstructuredGrid::SafeDownCast(grid);
1782   if(ugrid == NULL)
1783     return;
1784
1785   for(int index=0; index<vtkMedUtilities::GetNumberOfSubEntity(
1786       array->GetEntity().GeometryType); index++)
1787     {
1788     subentity.GeometryType = vtkMedUtilities::GetSubGeometry(
1789         array->GetEntity().GeometryType, index);
1790     vtkMedEntityArray* subarray=ugrid->GetEntityArray(subentity);
1791
1792     if(subarray==NULL)
1793       {
1794       vtkErrorMacro("DESC connectivity used, but sub types do not exist in file.");
1795       return;
1796       }
1797     subarray->LoadConnectivity();
1798     }
1799 }
1800
1801 vtkDataSet* vtkMedReader::CreateUnstructuredGridForPointSupport(
1802     vtkMedFamilyOnEntityOnProfile* foep)
1803 {
1804   vtkUnstructuredGrid* vtkgrid = vtkUnstructuredGrid::New();
1805   foep->ComputeIntersection(NULL);
1806   vtkMedFamilyOnEntity* foe = foep->GetFamilyOnEntity();
1807   vtkMedGrid* medgrid=foe->GetParentGrid();
1808   vtkMedUnstructuredGrid* medugrid=vtkMedUnstructuredGrid::SafeDownCast(
1809       medgrid);
1810   vtkMedCurvilinearGrid* medcgrid=vtkMedCurvilinearGrid::SafeDownCast(
1811       medgrid);
1812
1813   medgrid->LoadCoordinates();
1814
1815   vtkIdType npts=medgrid->GetNumberOfPoints();
1816
1817   bool shallowCopy= (medugrid != NULL || medcgrid!=NULL);
1818   if(medgrid->GetParentMesh()->GetSpaceDimension()!=3)
1819     {
1820     shallowCopy=false;
1821     }
1822   else
1823     {
1824     shallowCopy = foep->CanShallowCopyPointField(NULL);
1825     }
1826
1827   vtkDataArray* coords = NULL;
1828
1829   if(medugrid != NULL)
1830     coords = medugrid->GetCoordinates();
1831   if(medcgrid != NULL)
1832     coords = medcgrid->GetCoordinates();
1833
1834
1835   vtkIdType numberOfPoints;
1836   vtkPoints* points=vtkPoints::New(coords->GetDataType());
1837   vtkgrid->SetPoints(points);
1838   points->Delete();
1839
1840   vtkIdTypeArray* pointGlobalIds=vtkIdTypeArray::New();
1841   pointGlobalIds->SetName("MED_POINT_ID");
1842   pointGlobalIds->SetNumberOfComponents(1);
1843   vtkgrid->GetPointData()->SetGlobalIds(pointGlobalIds);
1844   pointGlobalIds->Delete();
1845
1846   if (shallowCopy)
1847     {
1848     vtkgrid->GetPoints()->SetData(coords);
1849     numberOfPoints=npts;
1850
1851     pointGlobalIds->SetNumberOfTuples(numberOfPoints);
1852     vtkIdType* ptr=pointGlobalIds->GetPointer(0);
1853     for(int pid=0; pid<numberOfPoints; pid++)
1854       ptr[pid]=pid+1;
1855     }
1856   if(!shallowCopy)
1857     {
1858     vtkIdType currentIndex=0;
1859
1860     for(vtkIdType index=0; index<medgrid->GetNumberOfPoints(); index++)
1861       {
1862       if (!foep->KeepPoint(index))
1863         {
1864         continue;
1865         }
1866
1867       double coord[3]={0.0, 0.0, 0.0};
1868       double * tuple=medgrid->GetCoordTuple(index);
1869       for(int dim=0; dim<medgrid->GetParentMesh()->GetSpaceDimension()&&dim<3; dim++)
1870         {
1871         coord[dim]=tuple[dim];
1872         }
1873       vtkgrid->GetPoints()->InsertPoint(currentIndex, coord);
1874       pointGlobalIds->InsertNextValue(index+1);
1875       currentIndex++;
1876       }
1877     vtkgrid->GetPoints()->Squeeze();
1878     pointGlobalIds->Squeeze();
1879     numberOfPoints=currentIndex;
1880     }
1881
1882   // now create the VTK_VERTEX cells
1883   for(vtkIdType id=0; id<numberOfPoints; id++)
1884     {
1885     vtkgrid->InsertNextCell(VTK_VERTEX, 1, &id);
1886     }
1887   vtkgrid->Squeeze();
1888
1889   // in this particular case, the global ids of the cells is the same as the global ids of the points.
1890   vtkgrid->GetCellData()->SetGlobalIds(vtkgrid->GetPointData()->GetGlobalIds());
1891
1892   return vtkgrid;
1893 }
1894
1895 vtkMedGrid* vtkMedReader::FindGridStep(vtkMedMesh* mesh)
1896 {
1897   if(this->Internal->RealAnimationMode == vtkMedReader::PhysicalTime)
1898     {
1899     vtkMedComputeStep cs;
1900     cs.TimeOrFrequency = this->Internal->UpdateTimeStep;
1901     return mesh->FindGridStep(cs, vtkMedReader::PhysicalTime);
1902     }
1903   else if(this->Internal->RealAnimationMode == vtkMedReader::Modes)
1904     {
1905     vtkMedComputeStep cs;
1906     cs.IterationIt = MED_NO_IT;
1907     cs.TimeIt = MED_NO_DT;
1908     cs.TimeOrFrequency = MED_NO_DT;
1909     return mesh->FindGridStep(cs, vtkMedReader::Modes);
1910     }
1911   else // Iterations
1912     {
1913     vtkMedComputeStep cs;
1914     // the time is set by choosing its index in the global
1915     // array giving the available times : this->TimeIndexForIterations
1916     cs.TimeOrFrequency = (med_int)this->AvailableTimes->GetValue(
1917         (vtkIdType)this->TimeIndexForIterations);
1918     // the iteration is asked by the pipeline
1919     cs.IterationIt = (med_int)this->Internal->UpdateTimeStep;
1920     return mesh->FindGridStep(cs, vtkMedReader::Iteration);
1921     }
1922   return NULL;
1923 }
1924
1925 void vtkMedReader::CreateMedSupports()
1926 {
1927   this->Internal->UsedSupports.clear();
1928
1929   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1930       meshfileit = this->Internal->MedFiles.begin();
1931   while(meshfileit != this->Internal->MedFiles.end())
1932     {
1933     vtkMedFile* file = meshfileit->second;
1934     meshfileit++;
1935
1936     for(int meshIndex=0; meshIndex<file->GetNumberOfMesh(); meshIndex++)
1937       {
1938       vtkMedMesh* mesh = file->GetMesh(meshIndex);
1939       vtkMedGrid* grid = this->FindGridStep(mesh);
1940       if(grid == NULL)
1941         continue;
1942
1943       for(int entityIndex=0; entityIndex<grid->GetNumberOfEntityArray();
1944         entityIndex++)
1945         {
1946         vtkMedEntityArray* array=grid->GetEntityArray(entityIndex);
1947         if(this->GetEntityStatus(array->GetEntity())==0)
1948           {
1949           continue;
1950           }
1951
1952         file->GetMedDriver()->LoadFamilyIds(array);
1953         for(int foeIndex=0; foeIndex<array->GetNumberOfFamilyOnEntity();
1954           foeIndex++)
1955           {
1956           vtkMedFamilyOnEntity* foe=array->GetFamilyOnEntity(foeIndex);
1957           vtkMedFamily* family=foe->GetFamily();
1958           if(this->GetFamilyStatus(mesh, family)==0)
1959             continue;
1960
1961           // now, I look over all non-point fields to see which profiles
1962           // have to be used on points.
1963           bool selectedSupport = false;
1964
1965           std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
1966               fieldfileit = this->Internal->MedFiles.begin();
1967           while(fieldfileit != this->Internal->MedFiles.end())
1968             {
1969             vtkMedFile* fieldfile = fieldfileit->second;
1970             fieldfileit++;
1971
1972             for(int fieldId=0; fieldId<fieldfile->GetNumberOfField(); fieldId++)
1973               {
1974               vtkMedField* field=fieldfile->GetField(fieldId);
1975
1976               if (!this->IsFieldSelected(field))
1977                 continue;
1978
1979               vtkMedListOfFieldSteps steps;
1980
1981               this->GatherFieldSteps(field, steps);
1982
1983               vtkMedListOfFieldSteps::iterator it=steps.begin();
1984               while(it!=steps.end())
1985                 {
1986                 vtkMedFieldStep *step = *it;
1987                 step->LoadInformation();
1988                 it++;
1989
1990                 for(int eid=0; eid<step->GetNumberOfFieldOverEntity(); eid++)
1991                   {
1992                   vtkMedFieldOverEntity* fieldOverEntity =
1993                       step->GetFieldOverEntity(eid);
1994
1995                   for(int pid = 0; pid < fieldOverEntity->GetNumberOfFieldOnProfile(); pid++)
1996                     {
1997                     vtkMedFieldOnProfile* fop =
1998                         fieldOverEntity->GetFieldOnProfile(pid);
1999                     vtkMedProfile* prof = fop->GetProfile();
2000                     vtkMedFamilyOnEntityOnProfile* foep =
2001                         foe->GetFamilyOnEntityOnProfile(prof);
2002                     if(foep != NULL)
2003                       {
2004                       this->Internal->UsedSupports.insert(foep);
2005                       selectedSupport = true;
2006                       }
2007                     }
2008                   }
2009                 }
2010               }
2011             }
2012           // If no field use this family on entity, I nevertheless create the
2013           // support, with an empty profile.
2014           if(!selectedSupport)
2015             {
2016             vtkMedFamilyOnEntityOnProfile* foep =
2017                 foe->GetFamilyOnEntityOnProfile((vtkMedProfile*)NULL);
2018             if(foep == NULL)
2019               {
2020               foep = vtkMedFamilyOnEntityOnProfile::New();
2021               foep->SetFamilyOnEntity(foe);
2022               foep->SetProfile(NULL);
2023               foe->AddFamilyOnEntityOnProfile(foep);
2024               foep->Delete();
2025               }
2026             this->Internal->UsedSupports.insert(foep);
2027             }
2028           }
2029         }
2030       }
2031   }
2032 }
2033
2034 bool vtkMedReader::BuildVTKSupport(
2035     vtkMedFamilyOnEntityOnProfile* foep,
2036     int doBuildSupport)
2037 {
2038
2039   vtkMedFamilyOnEntity* foe = foep->GetFamilyOnEntity();
2040
2041   int numProc = 1;
2042   vtkMultiProcessController* controller = vtkMultiProcessController::GetGlobalController();
2043   if (controller != NULL)
2044     {
2045     numProc = controller->GetNumberOfProcesses();
2046     }
2047
2048   if ((foep->GetValid() == 0) && numProc == 1)
2049     {
2050     return false;
2051     }
2052
2053   vtkMedGrid* grid=foe->GetParentGrid();
2054
2055   vtkMedEntityArray* array=foe->GetEntityArray();
2056   vtkMedMesh* mesh=grid->GetParentMesh();
2057   vtkSmartPointer<vtkStringArray> path = vtkSmartPointer<vtkStringArray>::New();
2058   string meshName=vtkMedUtilities::SimplifyName(mesh->GetName());
2059   path->InsertNextValue(meshName);
2060   std::string finalName;
2061
2062   if(foe->GetPointOrCell()==vtkMedUtilities::OnPoint)
2063     {
2064     path->InsertNextValue(vtkMedUtilities::OnPointName);
2065     finalName=vtkMedUtilities::SimplifyName(foe->GetFamily()->GetName());
2066     }
2067   else
2068     {
2069     path->InsertNextValue(vtkMedUtilities::OnCellName);
2070     path->InsertNextValue(vtkMedUtilities::SimplifyName(foe->GetFamily()->GetName()));
2071     finalName=vtkMedUtilities::EntityKey(array->GetEntity());
2072     }
2073
2074   if(foep->GetProfile() != NULL)
2075     {
2076     path->InsertNextValue(finalName);
2077     finalName = foep->GetProfile()->GetName();
2078     }
2079
2080   ostringstream progressBarTxt;
2081   for(int depth=0; depth<path->GetNumberOfValues(); depth++)
2082     {
2083     progressBarTxt<<path->GetValue(depth)<<" ";
2084     }
2085   progressBarTxt<<finalName;
2086   SetProgressText(progressBarTxt.str().c_str());
2087
2088   vtkDataSet* cachedDataSet = NULL;
2089   if(this->Internal->DataSetCache.find(foep)!=this->Internal->DataSetCache.end())
2090     {
2091     cachedDataSet = this->Internal->DataSetCache[foep];
2092     }
2093   else
2094     {
2095     vtkDataSet* dataset = NULL;
2096     if(doBuildSupport)
2097       {
2098       if(foe->GetPointOrCell()==vtkMedUtilities::OnPoint)
2099         {
2100         dataset = this->CreateUnstructuredGridForPointSupport(foep);
2101         }
2102       else
2103         {
2104         dataset = foep->GetFamilyOnEntity()->GetParentGrid()->
2105                   CreateVTKDataSet(foep);
2106         }
2107       }
2108
2109     if(dataset == NULL)
2110       {
2111       return false;
2112       }
2113
2114     this->Internal->DataSetCache[foep]=dataset;
2115     cachedDataSet = dataset;
2116     if(dataset != NULL)
2117       dataset->Delete();
2118   }
2119
2120   vtkMultiBlockDataSet* root=vtkMedUtilities::GetParent(this->GetOutput(), path);
2121   int nb=root->GetNumberOfBlocks();
2122
2123   if(cachedDataSet != NULL)
2124     {
2125     vtkDataSet* realDataSet=cachedDataSet->NewInstance();
2126     root->SetBlock(nb, realDataSet);
2127     realDataSet->Delete();
2128
2129     root->GetMetaData(nb)->Set(vtkCompositeDataSet::NAME(), finalName.c_str());
2130     realDataSet->ShallowCopy(cachedDataSet);
2131
2132     this->Internal->DataSetCache[foep]=cachedDataSet;
2133     this->Internal->CurrentDataSet[foep]=realDataSet;
2134
2135     path->InsertNextValue(finalName);
2136     path->SetName("BLOCK_NAME");
2137     realDataSet->GetFieldData()->AddArray(path);
2138     realDataSet->GetInformation()->Remove(vtkMedUtilities::BLOCK_NAME());
2139     for(int depth=0; depth<path->GetNumberOfValues(); depth++)
2140       {
2141       realDataSet->GetInformation()->Set(vtkMedUtilities::BLOCK_NAME(),
2142                                          path->GetValue(depth), depth);
2143       }
2144     }
2145 }
2146
2147 void vtkMedReader::MapFieldOnSupport(vtkMedFieldOnProfile* fop,
2148                                      vtkMedFamilyOnEntityOnProfile* foep,
2149                                      int doCreateField)
2150 {
2151   bool cached = false;
2152
2153   if(this->Internal->FieldCache.find(foep)
2154       !=this->Internal->FieldCache.end())
2155     {
2156     map<vtkMedFieldOnProfile*, VTKField>& fieldCache =
2157         this->Internal->FieldCache[foep];
2158     if(fieldCache.find(fop)!=fieldCache.end())
2159       {
2160       cached=true;
2161       }
2162     }
2163
2164   if(!cached)
2165     {
2166     this->CreateVTKFieldOnSupport(fop, foep, doCreateField);
2167     }
2168   this->SetVTKFieldOnSupport(fop, foep);
2169 }
2170
2171 void vtkMedReader::MapFieldsOnSupport(vtkMedFamilyOnEntityOnProfile* foep,
2172                                       int doCreateField)
2173 {
2174   // now loop over all fields to map it to the created grids
2175   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
2176       fieldfileit = this->Internal->MedFiles.begin();
2177   while(fieldfileit != this->Internal->MedFiles.end())
2178     {
2179     vtkMedFile* file = fieldfileit->second;
2180     fieldfileit++;
2181
2182     for(int fieldId=0; fieldId<file->GetNumberOfField(); fieldId++)
2183       {
2184       vtkMedField* field=file->GetField(fieldId);
2185
2186       if(strcmp(foep->GetFamilyOnEntity()->
2187                 GetParentGrid()->GetParentMesh()->GetName(),
2188                 field->GetMeshName()) != 0)
2189         continue;
2190
2191       if(strcmp(foep->GetFamilyOnEntity()->
2192                 GetParentGrid()->GetParentMesh()->GetName(),
2193                 field->GetMeshName()) != 0)
2194         continue;
2195
2196       if (!this->IsFieldSelected(field))
2197         continue;
2198
2199       vtkMedListOfFieldSteps steps;
2200
2201       this->GatherFieldSteps(field, steps);
2202
2203       vtkMedListOfFieldSteps::iterator it=steps.begin();
2204       while(it!=steps.end())
2205         {
2206         vtkMedFieldStep *step = *it;
2207         step->LoadInformation();
2208         it++;
2209
2210         for(int eid=0; eid<step->GetNumberOfFieldOverEntity(); eid++)
2211           {
2212           vtkMedFieldOverEntity* fieldOverEntity = step->GetFieldOverEntity(eid);
2213           for(int pid = 0; pid < fieldOverEntity->GetNumberOfFieldOnProfile(); pid++)
2214             {
2215             vtkMedFieldOnProfile* fop = fieldOverEntity->GetFieldOnProfile(pid);
2216             if(foep->CanMapField(fop))
2217               {
2218               this->MapFieldOnSupport(fop, foep, doCreateField);
2219               }
2220             }
2221           }
2222         }
2223       }
2224     }
2225 }
2226
2227 void vtkMedReader::GatherFieldSteps(vtkMedField* field,
2228                                     vtkMedListOfFieldSteps& steps)
2229 {
2230   if(this->Internal->RealAnimationMode == vtkMedReader::PhysicalTime)
2231     {
2232     vtkMedComputeStep cs;
2233     cs.TimeOrFrequency = this->Internal->UpdateTimeStep;
2234     vtkMedFieldStep* fs =
2235         field->FindFieldStep(cs, vtkMedReader::PhysicalTime);
2236     if(fs != NULL)
2237       steps.push_back(fs);
2238     }
2239   else if(this->Internal->RealAnimationMode == vtkMedReader::Iteration)
2240     {
2241     vtkMedComputeStep cs;
2242     cs.IterationIt = (med_int)this->Internal->UpdateTimeStep;
2243     cs.TimeOrFrequency = (med_int)this->AvailableTimes->GetValue(
2244         (vtkIdType)this->TimeIndexForIterations);
2245     vtkMedFieldStep* fs =
2246         field->FindFieldStep(cs, vtkMedReader::Iteration);
2247     if(fs != NULL)
2248       steps.push_back(fs);
2249     }
2250   else // modes
2251     {
2252     for(int modeid = 0; modeid < this->Frequencies->GetNumberOfKey(); modeid++)
2253       {
2254       if(this->Frequencies->GetKeyStatus(
2255           this->Frequencies->GetKey(modeid)) == 0)
2256         {
2257         continue;
2258         }
2259
2260       vtkMedComputeStep cs;
2261       cs.TimeOrFrequency = this->AvailableTimes->GetValue(modeid);
2262       vtkMedFieldStep* fs =
2263           field->FindFieldStep(cs, vtkMedReader::PhysicalTime);
2264       if(fs != NULL)
2265         steps.push_back(fs);
2266       }
2267     }
2268 }
2269
2270 void vtkMedReader::SetVTKFieldOnSupport(vtkMedFieldOnProfile* fop,
2271     vtkMedFamilyOnEntityOnProfile* foep)
2272 {
2273   vtkMedFamilyOnEntity* foe = foep->GetFamilyOnEntity();
2274   vtkMedFieldOverEntity* fieldOverEntity = fop->GetParentFieldOverEntity();
2275   vtkMedFieldStep* step = fieldOverEntity->GetParentStep();
2276   vtkMedField* field = step->GetParentField();
2277
2278   const vtkMedComputeStep& cs = step->GetComputeStep();
2279
2280   vtkDataSet* ds=this->Internal->CurrentDataSet[foep];
2281   if (ds == NULL)
2282     {
2283   // ds == NULL could arrive is some cases when working in parallel
2284   vtkWarningMacro( "--- vtkMedReader::SetVTKFieldOnSupport: ds is NULL !!");
2285   return;
2286     }
2287
2288   VTKField& vtkfield=this->Internal->FieldCache[foep][fop];
2289
2290   std::string name=vtkMedUtilities::SimplifyName(field->GetName());
2291   std::string vectorname = name+"_Vector";
2292   if(this->AnimationMode==Modes)
2293     {
2294     double freq = cs.TimeOrFrequency;
2295     int index = this->GetFrequencyIndex(freq);
2296     name += string(" ") + vtkMedUtilities::GetModeKey(index, freq,
2297       this->AvailableTimes->GetNumberOfTuples()-1);
2298     vectorname += string(" ") + vtkMedUtilities::GetModeKey(index, freq,
2299       this->AvailableTimes->GetNumberOfTuples()-1);
2300     }
2301
2302   vtkfield.DataArray->SetName(name.c_str());
2303   vtkfield.DataArray->Squeeze();
2304   if(vtkfield.Vectors != NULL)
2305     {
2306     vtkfield.Vectors->SetName(vectorname.c_str());
2307     vtkfield.Vectors->Squeeze();
2308     }
2309   if(vtkfield.QuadratureIndexArray!=NULL)
2310     {
2311     vtkfield.QuadratureIndexArray->Squeeze();
2312     }
2313
2314   if(foe->GetPointOrCell()==vtkMedUtilities::OnCell)
2315     {
2316     if(field->GetFieldType()==vtkMedField::PointField)
2317       {
2318       if(vtkfield.DataArray->GetNumberOfTuples()!=ds->GetNumberOfPoints())
2319         {
2320           vtkDebugMacro("the data array " << vtkfield.DataArray->GetName()
2321                       << " do not have the good number of tuples");
2322         return;
2323         }
2324       ds->GetPointData()->AddArray(vtkfield.DataArray);
2325       if(vtkfield.Vectors != NULL && this->GenerateVectors)
2326         {
2327         ds->GetPointData()->AddArray(vtkfield.Vectors);
2328         ds->GetPointData()->SetActiveVectors(vtkfield.Vectors->GetName());
2329         }
2330       switch (vtkfield.DataArray->GetNumberOfComponents())
2331         {
2332         case 1:
2333         ds->GetPointData()->SetActiveScalars(vtkfield.DataArray->GetName());
2334         break;
2335         case 3:
2336         ds->GetPointData()->SetActiveVectors(vtkfield.DataArray->GetName());
2337         break;
2338         }
2339       // if the data set is only composed of VTK_VERTEX cells,
2340       // and no field called with the same name exist on cells,
2341       // map this field to cells too
2342       if(foe->GetVertexOnly()==1&&ds->GetCellData()->GetArray(
2343               vtkfield.DataArray->GetName())==NULL)
2344         {
2345         ds->GetCellData()->AddArray(vtkfield.DataArray);
2346         if(vtkfield.Vectors != NULL && this->GenerateVectors)
2347           {
2348           ds->GetCellData()->AddArray(vtkfield.Vectors);
2349           ds->GetCellData()->SetActiveVectors(vtkfield.Vectors->GetName());
2350           }
2351         switch (vtkfield.DataArray->GetNumberOfComponents())
2352           {
2353           case 1:
2354           ds->GetCellData()->SetActiveScalars(vtkfield.DataArray->GetName());
2355           break;
2356           case 3:
2357           ds->GetCellData()->SetActiveVectors(vtkfield.DataArray->GetName());
2358           break;
2359           }
2360         }
2361       }
2362     if(field->GetFieldType()==vtkMedField::CellField)
2363       {
2364       if((this->Internal->NumberOfPieces == 1) && vtkfield.DataArray->GetNumberOfTuples()!=ds->GetNumberOfCells()  )
2365         {
2366         vtkDebugMacro("the data array " << vtkfield.DataArray->GetName()
2367                       << " do not have the good number of tuples"
2368                       << " ncell=" << ds->GetNumberOfCells()
2369                       << " ntuple=" << vtkfield.DataArray->GetNumberOfTuples());
2370         return;
2371         }
2372       // In case we are in parallel and our process does not contain the data
2373       if(ds->GetNumberOfCells()==0)
2374         {
2375         return;
2376         }
2377       ds->GetCellData()->AddArray(vtkfield.DataArray);
2378
2379       if(vtkfield.Vectors != NULL && this->GenerateVectors)
2380         {
2381         ds->GetCellData()->AddArray(vtkfield.Vectors);
2382         ds->GetCellData()->SetActiveVectors(vtkfield.Vectors->GetName());
2383         }
2384       switch (vtkfield.DataArray->GetNumberOfComponents())
2385         {
2386         case 1:
2387         ds->GetCellData()->SetActiveScalars(vtkfield.DataArray->GetName());
2388         break;
2389         case 3:
2390         ds->GetCellData()->SetActiveVectors(vtkfield.DataArray->GetName());
2391         break;
2392         }
2393       // if the data set is only composed of VTK_VERTEX cells,
2394       // and no field called with the same name exist on points,
2395       // map this field to points too
2396       if(foe->GetVertexOnly()==1 && ds->GetPointData()->GetArray(
2397               vtkfield.DataArray->GetName())==NULL)
2398         {
2399         ds->GetPointData()->AddArray(vtkfield.DataArray);
2400         if(vtkfield.Vectors != NULL && this->GenerateVectors)
2401           {
2402           ds->GetPointData()->AddArray(vtkfield.Vectors);
2403           ds->GetPointData()->SetActiveVectors(vtkfield.Vectors->GetName());
2404           }
2405         switch (vtkfield.DataArray->GetNumberOfComponents())
2406           {
2407           case 1:
2408           ds->GetPointData()->SetActiveScalars(vtkfield.DataArray->GetName());
2409           break;
2410           case 3:
2411           ds->GetPointData()->SetActiveVectors(vtkfield.DataArray->GetName());
2412           break;
2413           }
2414         }
2415       }
2416     if(field->GetFieldType()==vtkMedField::QuadratureField ||
2417        field->GetFieldType()==vtkMedField::ElnoField )
2418       {
2419       vtkIdType ncells=ds->GetNumberOfCells();
2420       vtkIdType nid=vtkfield.QuadratureIndexArray->GetNumberOfTuples();
2421       vtkIdType nda=vtkfield.DataArray->GetNumberOfTuples();
2422       if((nid!=ncells) && (this->Internal->NumberOfPieces == 1))
2423         {
2424         vtkDebugMacro(
2425             "There should be as many quadrature index values as there are cells");
2426         return;
2427         }
2428       else
2429         {
2430       if (ncells == 0)
2431       {
2432         vtkfield.DataArray->SetNumberOfTuples( 0 );
2433         vtkfield.DataArray->Squeeze();
2434       }
2435       if (nid>ncells)  // PROBABLY NOT NECESSARY
2436       {
2437         vtkfield.QuadratureIndexArray->SetNumberOfTuples(ncells);
2438         int nquad=fop->GetNumberOfIntegrationPoint();
2439         vtkfield.DataArray->SetNumberOfTuples( nquad * ds->GetNumberOfCells() );
2440         vtkfield.DataArray->Squeeze();
2441       }
2442         ds->GetFieldData()->AddArray(vtkfield.DataArray);
2443         ds->GetCellData()->AddArray(vtkfield.QuadratureIndexArray);
2444
2445         nid=vtkfield.QuadratureIndexArray->GetNumberOfTuples();
2446         nda=vtkfield.DataArray->GetNumberOfTuples();
2447
2448         if(vtkfield.Vectors != NULL && this->GenerateVectors)
2449           {
2450           ds->GetFieldData()->AddArray(vtkfield.Vectors);
2451           }
2452         }
2453       }
2454     }//support OnCell
2455   else
2456     {//support OnPoint
2457     if(vtkfield.DataArray->GetNumberOfTuples()!=ds->GetNumberOfPoints())
2458       {
2459       vtkDebugMacro("the data array " << vtkfield.DataArray->GetName() << " do not have the good number of tuples");
2460       return;
2461       }
2462     ds->GetPointData()->AddArray(vtkfield.DataArray);
2463     if(vtkfield.Vectors != NULL && this->GenerateVectors)
2464       {
2465       ds->GetPointData()->AddArray(vtkfield.Vectors);
2466       ds->GetPointData()->SetActiveVectors(vtkfield.Vectors->GetName());
2467       }
2468     switch (vtkfield.DataArray->GetNumberOfComponents())
2469       {
2470       case 1:
2471       ds->GetPointData()->SetActiveScalars(vtkfield.DataArray->GetName());
2472       break;
2473       case 3:
2474       ds->GetPointData()->SetActiveVectors(vtkfield.DataArray->GetName());
2475       break;
2476       }
2477     // all the VTK_VERTEX created cells have the same order than the points,
2478     // I can safely map the point array to the cells in this particular case.
2479     ds->GetCellData()->AddArray(vtkfield.DataArray);
2480     if(vtkfield.Vectors != NULL && this->GenerateVectors)
2481       {
2482       ds->GetCellData()->AddArray(vtkfield.Vectors);
2483       ds->GetCellData()->SetActiveVectors(vtkfield.Vectors->GetName());
2484       }
2485     switch (vtkfield.DataArray->GetNumberOfComponents())
2486       {
2487       case 1:
2488       ds->GetCellData()->SetActiveScalars(vtkfield.DataArray->GetName());
2489       break;
2490       case 3:
2491       ds->GetCellData()->SetActiveVectors(vtkfield.DataArray->GetName());
2492       break;
2493       }
2494     }
2495 }
2496
2497 void vtkMedReader::InitializeQuadratureOffsets(vtkMedFieldOnProfile* fop,
2498     vtkMedFamilyOnEntityOnProfile* foep)
2499 {
2500   vtkMedFamilyOnEntity* foe = foep->GetFamilyOnEntity();
2501
2502   vtkMedFieldOverEntity* fieldOverEntity = fop->GetParentFieldOverEntity();
2503   vtkMedFieldStep* step = fieldOverEntity->GetParentStep();
2504   vtkMedField *field= step->GetParentField();
2505
2506   // then I compute the quadrature key if needed, and try to find the quadrature offsets.
2507   if(this->Internal->QuadratureOffsetCache.find(foep)
2508       ==this->Internal->QuadratureOffsetCache.end())
2509     this->Internal->QuadratureOffsetCache[foep]=map<LocalizationKey,
2510         vtkSmartPointer<vtkIdTypeArray> > ();
2511
2512   map<LocalizationKey, vtkSmartPointer<vtkIdTypeArray> >& quadOffsets=
2513       this->Internal->QuadratureOffsetCache[foep];
2514
2515   LocalizationKey quadkey=this->GetLocalizationKey(fop);
2516
2517   if(quadOffsets.find(quadkey)!=quadOffsets.end())
2518     {// the quadrature offset array has already been created
2519     return;
2520     }
2521
2522   vtkIdTypeArray* qoffsets=vtkIdTypeArray::New();
2523   quadOffsets[quadkey]=qoffsets;
2524   qoffsets->Delete();
2525
2526   ostringstream sstr;
2527   if(field->GetFieldType() == vtkMedField::ElnoField)
2528     {
2529     qoffsets->GetInformation()->Set(vtkMedUtilities::ELNO(), 1);
2530     sstr<<"ELNO";
2531     }
2532   else if(field->GetFieldType() == vtkMedField::QuadratureField)
2533     {
2534     qoffsets->GetInformation()->Set(vtkMedUtilities::ELGA(), 1);
2535     sstr<<"ELGA";
2536     }
2537   else
2538     {
2539     sstr<<"QuadraturePointOffset";
2540     }
2541
2542   qoffsets->SetName(sstr.str().c_str());
2543
2544   vtkSmartPointer<vtkMedLocalization> loc=
2545       this->GetLocalization(fop->GetLocalizationName());
2546
2547   if(loc == NULL)
2548     {
2549     if(fop->GetNumberOfIntegrationPoint()==1)
2550       {// cell-centered fields can be seen as quadrature fields with 1
2551       // quadrature point at the center
2552       vtkMedLocalization* center=vtkMedLocalization::New();
2553       loc=center;
2554       center->Delete();
2555       center->BuildCenter(fieldOverEntity->GetEntity().GeometryType);
2556       }
2557     else if(loc == NULL && field->GetFieldType() == vtkMedField::ElnoField)
2558       {// ELNO fields have no vtkMedLocalization,
2559       // I need to create a dummy one
2560       vtkMedLocalization* elnodef=vtkMedLocalization::New();
2561       loc=elnodef;
2562       elnodef->Delete();
2563       elnodef->BuildELNO(fieldOverEntity->GetEntity().GeometryType);
2564       }
2565     else
2566       {
2567       vtkErrorMacro("Cannot find localization of quadrature field "
2568                     << field->GetName());
2569       }
2570     }
2571   this->AddQuadratureSchemeDefinition(qoffsets->GetInformation(), loc);
2572 }
2573
2574 void vtkMedReader::CreateVTKFieldOnSupport(vtkMedFieldOnProfile* fop,
2575     vtkMedFamilyOnEntityOnProfile* foep, int doCreateField)
2576 {
2577   vtkMedFamilyOnEntity* foe = foep->GetFamilyOnEntity();
2578   vtkMedFieldOverEntity* fieldOverEntity = fop->GetParentFieldOverEntity();
2579   vtkMedFieldStep* step = fieldOverEntity->GetParentStep();
2580   vtkMedField* field = step->GetParentField();
2581
2582   if(vtkMedUnstructuredGrid::SafeDownCast(
2583       foep->GetFamilyOnEntity()->GetParentGrid()) != NULL)
2584     {
2585     fop->Load(MED_COMPACT_STMODE);
2586     }
2587   else
2588     {
2589     fop->Load(MED_GLOBAL_STMODE);
2590     }
2591
2592   vtkMedIntArray* profids=NULL;
2593
2594   vtkMedProfile* profile=fop->GetProfile();
2595
2596   if(profile)
2597     {
2598     profile->Load();
2599     profids=profile->GetIds();
2600     }//has profile
2601
2602   VTKField& vtkfield=this->Internal->FieldCache[foep][fop];
2603
2604   bool shallowok = true;
2605
2606   // for structured grid, the values are loaded globally, and cells which are
2607   // not on the profile or not on the family are blanked out.
2608   // shallow copy is therefore always possible
2609   if(vtkMedUnstructuredGrid::SafeDownCast(
2610       foep->GetFamilyOnEntity()->GetParentGrid()) != NULL)
2611     {
2612     shallowok = foep->CanShallowCopy(fop);
2613     }
2614
2615   if(shallowok)
2616     {
2617     vtkfield.DataArray = fop->GetData();
2618     }
2619   else
2620     {
2621     vtkDataArray* data=vtkMedUtilities::NewArray(field->GetDataType());
2622     vtkfield.DataArray=data;
2623     data->Delete();
2624     vtkfield.DataArray->SetNumberOfComponents(field->GetNumberOfComponent());
2625     }
2626
2627   for(int comp=0; comp<field->GetNumberOfComponent(); comp++)
2628     {
2629     vtkfield.DataArray->SetComponentName(comp, vtkMedUtilities::SimplifyName(
2630         field->GetComponentName()->GetValue(comp)).c_str());
2631     }
2632
2633   bool createOffsets=false;
2634   if(field->GetFieldType()==vtkMedField::QuadratureField ||
2635      field->GetFieldType()==vtkMedField::ElnoField )
2636     {
2637     this->InitializeQuadratureOffsets(fop, foep);
2638
2639     LocalizationKey quadKey = this->GetLocalizationKey(fop);
2640     vtkfield.QuadratureIndexArray
2641         =this->Internal->QuadratureOffsetCache[foep][quadKey];
2642     vtkDataSet* ds = this->Internal->CurrentDataSet[foep];
2643
2644     vtkfield.DataArray->GetInformation()->Set(
2645         vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME(),
2646         vtkfield.QuadratureIndexArray->GetName());
2647     vtkfield.QuadratureIndexArray->GetInformation()->Set(
2648         vtkAbstractArray::GUI_HIDE(), 1);
2649
2650     if(vtkfield.QuadratureIndexArray->GetNumberOfTuples()
2651         !=ds->GetNumberOfCells())
2652       {
2653       vtkfield.QuadratureIndexArray->SetNumberOfTuples(0);
2654       createOffsets=true;
2655       }
2656     }
2657
2658   if(!doCreateField)
2659     return;
2660
2661   if(shallowok)
2662     {
2663     // build the quadrature offset array if needed
2664     if(createOffsets)
2665       {
2666       vtkIdType noffsets;
2667       int nquad=fop->GetNumberOfIntegrationPoint();
2668       noffsets=fop->GetData()->GetNumberOfTuples()/nquad;
2669
2670       vtkDataSet* ds=this->Internal->CurrentDataSet[foep];
2671       if(noffsets!=ds->GetNumberOfCells())
2672         {
2673         vtkErrorMacro(
2674             "number of quadrature offsets (" << noffsets << ") must "
2675             << "match the number of cells (" << ds->GetNumberOfCells() << ")!");
2676         }
2677
2678       vtkfield.QuadratureIndexArray->SetNumberOfTuples(noffsets);
2679       for(vtkIdType id=0; id<noffsets; id++)
2680         {
2681         vtkfield.QuadratureIndexArray->SetValue(id, nquad*id);
2682         }
2683       }
2684
2685     }
2686   else if(foe->GetPointOrCell() == vtkMedUtilities::OnCell
2687      && field->GetFieldType() != vtkMedField::PointField)
2688     {
2689     // Cell-centered field on cell support
2690     int nquad = 1;
2691     if (field->GetFieldType()==vtkMedField::QuadratureField ||
2692         field->GetFieldType()==vtkMedField::ElnoField)
2693       {
2694       nquad=fop->GetNumberOfIntegrationPoint();
2695       }
2696     vtkMedIntArray* profids=NULL;
2697     if (profile)
2698       {
2699       profids=profile->GetIds();
2700       }
2701     vtkIdType maxIndex=fop->GetData()->GetNumberOfTuples();
2702     maxIndex/=nquad;
2703     vtkIdType quadIndex = 0;
2704
2705     for (vtkIdType id = 0; id<maxIndex; id++)
2706       {
2707       vtkIdType realIndex = (profids!=NULL ? profids->GetValue(id)-1 : id);
2708       if (!foep->KeepCell(realIndex))
2709         continue;
2710
2711       if (field->GetFieldType()==vtkMedField::QuadratureField ||
2712           field->GetFieldType()==vtkMedField::ElnoField)
2713         {
2714         for (int q = 0; q<nquad; q++)
2715           {
2716           vtkfield.DataArray->InsertNextTuple(nquad*realIndex+q,
2717               fop->GetData());
2718           }
2719         if (createOffsets)
2720           {
2721           vtkfield.QuadratureIndexArray->InsertNextValue(quadIndex);
2722           quadIndex+=nquad;
2723           }
2724         }
2725       else
2726         {
2727         vtkfield.DataArray->InsertNextTuple(id, fop->GetData());
2728         }
2729       }//copy all tuples
2730     vtkfield.DataArray->Squeeze();
2731     vtkDataSet* ds = this->Internal->CurrentDataSet[foep];
2732     }
2733   else if(foe->GetPointOrCell() == vtkMedUtilities::OnCell
2734      && field->GetFieldType() == vtkMedField::PointField)
2735     {// point field on cell support
2736     vtkMedIntArray* profids=NULL;
2737
2738     vtkMedProfile* profile=fop->GetProfile();
2739
2740     if(profile)
2741       {
2742       profile->Load();
2743
2744       profids=profile->GetIds();
2745       }//has profile
2746
2747     vtkIdType maxId=fop->GetData()->GetNumberOfTuples();
2748
2749     for(vtkIdType id=0; id<maxId; id++)
2750       {
2751       // if I have a profile, then I should insert the value at the position given in the profile.
2752       vtkIdType destIndex;
2753       if(profids!=NULL)
2754         {
2755         destIndex=profids->GetValue(id)-1; // -1 because med indices start at 1
2756         }
2757       else
2758         {
2759         destIndex=id;
2760         }
2761
2762       if(!foep->KeepPoint(destIndex))
2763         continue;
2764       // if I use the med2VTKIndex, then the index to insert
2765       // this value is given by the map.
2766       destIndex = foep->GetVTKPointIndex(destIndex);
2767       vtkfield.DataArray->InsertTuple(destIndex, id, fop->GetData());
2768       }
2769     vtkfield.DataArray->Squeeze();
2770     }// point field on cell support
2771   else if(foe->GetPointOrCell() == vtkMedUtilities::OnPoint &&
2772      field->GetFieldType() == vtkMedField::PointField)
2773     {//support OnPoint
2774
2775     vtkIdType maxId = fop->GetData()->GetNumberOfTuples();
2776
2777     for(vtkIdType id=0; id<maxId; id++)
2778       {
2779       vtkIdType realIndex=id;
2780       if(profids!=NULL)
2781         {
2782         realIndex=profids->GetValue(id)-1; // -1 because med indices start at 1
2783         }
2784
2785       if(!foep->KeepPoint(realIndex))
2786         continue;
2787
2788       vtkfield.DataArray->InsertNextTuple(fop->GetData()->GetTuple(realIndex));
2789       }
2790     vtkfield.DataArray->Squeeze();
2791     }// support on point
2792
2793   // now generate the vector field if asked for
2794   if(this->GenerateVectors)
2795     {
2796     int ncomp = vtkfield.DataArray->GetNumberOfComponents();
2797     if(ncomp > 1 && ncomp != 3)
2798       {
2799       vtkDataArray* vectors = vtkfield.DataArray->NewInstance();
2800       vectors->SetNumberOfComponents(3);
2801       vectors->SetNumberOfTuples(vtkfield.DataArray->GetNumberOfTuples());
2802       vtkfield.Vectors = vectors;
2803       vectors->Delete();
2804
2805       vectors->CopyInformation(vtkfield.DataArray->GetInformation());
2806
2807       if(ncomp < 3)
2808         vectors->SetComponentName(2, "Z");
2809       else
2810         vectors->SetComponentName(2, vtkfield.DataArray->GetComponentName(2));
2811
2812       vectors->SetComponentName(1, vtkfield.DataArray->GetComponentName(1));
2813       vectors->SetComponentName(0, vtkfield.DataArray->GetComponentName(0));
2814
2815       int tuplesize = (ncomp > 3? ncomp: 3);
2816       double* tuple = new double[tuplesize];
2817       for(int tid=0; tid<tuplesize; tid++)
2818         {
2819         tuple[tid] = 0.0;
2820         }
2821
2822       for(vtkIdType id=0; id < vtkfield.DataArray->GetNumberOfTuples(); id++)
2823         {
2824         vtkfield.DataArray->GetTuple(id, tuple);
2825         vectors->SetTuple(id, tuple);
2826         }
2827       }
2828     }
2829 }
2830
2831 int vtkMedReader::HasMeshAnyCellSelectedFamily(vtkMedMesh* mesh)
2832 {
2833   int nfam = mesh->GetNumberOfCellFamily();
2834   for (int famid = 0; famid<nfam; famid++)
2835     {
2836     vtkMedFamily* fam = mesh->GetFamily(famid);
2837     if (fam->GetPointOrCell()!=vtkMedUtilities::OnCell||!this->GetFamilyStatus(
2838         mesh, fam))
2839       continue;
2840     return true;
2841     }
2842   return false;
2843 }
2844
2845 int vtkMedReader::HasMeshAnyPointSelectedFamily(vtkMedMesh* mesh)
2846 {
2847   int nfam = mesh->GetNumberOfCellFamily();
2848   for (int famid = 0; famid<nfam; famid++)
2849     {
2850     vtkMedFamily* fam = mesh->GetFamily(famid);
2851     if (fam->GetPointOrCell()!=vtkMedUtilities::OnPoint
2852         ||!this->GetFamilyStatus(mesh, fam))
2853       continue;
2854     return true;
2855     }
2856   return false;
2857 }
2858
2859 void vtkMedReader::BuildSIL(vtkMutableDirectedGraph* sil)
2860 {
2861   if(sil==NULL)
2862     return;
2863
2864   sil->Initialize();
2865
2866   vtkSmartPointer<vtkVariantArray> childEdge=
2867       vtkSmartPointer<vtkVariantArray>::New();
2868   childEdge->InsertNextValue(0);
2869
2870   vtkSmartPointer<vtkVariantArray> crossEdge=
2871       vtkSmartPointer<vtkVariantArray>::New();
2872   crossEdge->InsertNextValue(1);
2873
2874   // CrossEdge is an edge linking hierarchies.
2875   vtkUnsignedCharArray* crossEdgesArray=vtkUnsignedCharArray::New();
2876   crossEdgesArray->SetName("CrossEdges");
2877   sil->GetEdgeData()->AddArray(crossEdgesArray);
2878   crossEdgesArray->Delete();
2879   /*vtk*/std::deque</*vtk*/std::string> names;
2880
2881   // Now build the hierarchy.
2882   vtkIdType rootId=sil->AddVertex();
2883   names.push_back("SIL");
2884
2885   // Add a global entry to encode global names for the families
2886   vtkIdType globalFamilyRoot=sil->AddChild(rootId, childEdge);
2887   names.push_back("Families");
2888
2889   // Add a global entry to encode global names for the families
2890   vtkIdType globalGroupRoot=sil->AddChild(rootId, childEdge);
2891   names.push_back("Groups");
2892
2893   // Add the groups subtree
2894   vtkIdType groupsRoot=sil->AddChild(rootId, childEdge);
2895   names.push_back("GroupTree");
2896
2897   // Add a global entry to encode names for the cell types
2898   vtkIdType globalEntityRoot=sil->AddChild(rootId, childEdge);
2899   names.push_back("Entity");
2900
2901   // Add the cell types subtree
2902   vtkIdType entityTypesRoot=sil->AddChild(rootId, childEdge);
2903   names.push_back("EntityTree");
2904
2905   // this is the map that keep added cell types
2906   map<vtkMedEntity, vtkIdType> entityMap;
2907   map<med_entity_type, vtkIdType> typeMap;
2908
2909   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
2910       meshfileit = this->Internal->MedFiles.begin();
2911   while(meshfileit != this->Internal->MedFiles.end())
2912     {
2913     vtkMedFile* file = meshfileit->second;
2914     meshfileit++;
2915
2916     int numMeshes=file->GetNumberOfMesh();
2917     for(int meshIndex=0; meshIndex<numMeshes; meshIndex++)
2918       {
2919       vtkMedMesh* mesh = file->GetMesh(meshIndex);
2920       vtkMedGrid* grid = this->FindGridStep(mesh);
2921
2922       if(grid == NULL)
2923         continue;
2924
2925       // add all entities
2926       for(int entityIndex=0; entityIndex<grid->GetNumberOfEntityArray(); entityIndex++)
2927         {
2928         vtkMedEntityArray* array=grid->GetEntityArray(entityIndex);
2929         vtkMedEntity entity = array->GetEntity();
2930
2931         if(entityMap.find(entity)==entityMap.end())
2932           {
2933           vtkIdType entityGlobalId=sil->AddChild(globalEntityRoot, childEdge);
2934           names.push_back(vtkMedUtilities::EntityKey(entity));
2935
2936           vtkIdType typeId;
2937           if(typeMap.find(entity.EntityType)==typeMap.end())
2938             {
2939             typeId=sil->AddChild(entityTypesRoot, childEdge);
2940             names.push_back(vtkMedUtilities::EntityName(entity.EntityType));
2941             typeMap[entity.EntityType]=typeId;
2942             }
2943           else
2944             {
2945             typeId=typeMap[entity.EntityType];
2946             }
2947           vtkIdType entityId=sil->AddChild(typeId, childEdge);
2948           names.push_back(entity.GeometryName);
2949
2950           sil->AddEdge(entityId, entityGlobalId, crossEdge);
2951
2952           entityMap[entity]=entityId;
2953           }
2954         }
2955
2956       vtkIdType meshGroup = sil->AddChild(groupsRoot, childEdge);
2957       names.push_back(vtkMedUtilities::SimplifyName(mesh->GetName()));
2958
2959       // add the two OnPoint and OnCell entries, for groups and for families
2960       vtkIdType meshCellGroups = sil->AddChild(meshGroup, childEdge);
2961       names.push_back(vtkMedUtilities::OnCellName);
2962
2963       vtkIdType meshPointGroups = sil->AddChild(meshGroup, childEdge);
2964       names.push_back(vtkMedUtilities::OnPointName);
2965
2966       // this maps will keep all added groups on nodes/cells of this mesh
2967       map<string, vtkIdType> nodeGroupMap;
2968       map<string, vtkIdType> cellGroupMap;
2969
2970       // add all families
2971       for(int famIndex=0; famIndex<mesh->GetNumberOfFamily(); famIndex++)
2972         {
2973         vtkMedFamily* family=mesh->GetFamily(famIndex);
2974
2975         vtkIdType globalFamilyId = sil->AddChild(globalFamilyRoot, childEdge);
2976         names.push_back(vtkMedUtilities::FamilyKey(mesh->GetName(),
2977                                                    family->GetPointOrCell(),
2978                                                    family->GetName()));
2979
2980         // family with Id 0 is both on cell and on points, so add it to both
2981         map<string, vtkIdType> & groupMap=(family->GetPointOrCell()
2982             ==vtkMedUtilities::OnPoint? nodeGroupMap: cellGroupMap);
2983
2984         vtkIdType groupRootId =
2985             (family->GetPointOrCell()==vtkMedUtilities::OnPoint?
2986              meshPointGroups : meshCellGroups);
2987
2988         // add all the groups of this family
2989         for(vtkIdType groupIndex=0; groupIndex<family->GetNumberOfGroup();
2990           groupIndex++)
2991           {
2992           vtkMedGroup* group=family->GetGroup(groupIndex);
2993
2994           vtkIdType familyGroupId = sil->AddChild(globalFamilyId, childEdge);
2995           names.push_back(vtkMedUtilities::FamilyKey(
2996               mesh->GetName(), family->GetPointOrCell(),
2997               family->GetName()));
2998
2999           vtkIdType groupGlobalId;
3000           if(groupMap.find(group->GetName())==groupMap.end())
3001             {
3002             vtkIdType groupLocalId;
3003             groupLocalId=sil->AddChild(groupRootId, childEdge);
3004             names.push_back(vtkMedUtilities::SimplifyName(group->GetName()));
3005
3006             groupGlobalId=sil->AddChild(globalGroupRoot, childEdge);
3007             names.push_back(vtkMedUtilities::GroupKey(
3008                 mesh->GetName(), family->GetPointOrCell(),
3009                 group->GetName()));
3010             groupMap[group->GetName()]=groupGlobalId;
3011
3012             sil->AddEdge(groupLocalId, groupGlobalId, crossEdge);
3013             }
3014           vtkIdType groupId = groupMap[group->GetName()];
3015           sil->AddEdge(familyGroupId, groupId, childEdge);
3016
3017           }//groupIndex
3018         }//famIndex
3019       }//meshIndex
3020     }// file iterator
3021
3022   // This array is used to assign names to nodes.
3023   vtkStringArray* namesArray=vtkStringArray::New();
3024   namesArray->SetName("Names");
3025   namesArray->SetNumberOfTuples(sil->GetNumberOfVertices());
3026   sil->GetVertexData()->AddArray(namesArray);
3027   namesArray->Delete();
3028   /*vtk*/std::deque</*vtk*/std::string>::iterator iter;
3029   vtkIdType cc;
3030   for(cc=0, iter=names.begin(); iter!=names.end(); ++iter, ++cc)
3031     {
3032     namesArray->SetValue(cc, (*iter).c_str());
3033     }
3034 }
3035
3036 void vtkMedReader::ClearMedSupports()
3037 {
3038   this->Internal->DataSetCache.clear();
3039   //this->Internal->Med2VTKPointIndex.clear();
3040
3041   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
3042       meshfileit = this->Internal->MedFiles.begin();
3043   while(meshfileit != this->Internal->MedFiles.end())
3044     {
3045     vtkMedFile* file = meshfileit->second;
3046     meshfileit++;
3047     int numMeshes=file->GetNumberOfMesh();
3048     for(int meshIndex=0; meshIndex<numMeshes; meshIndex++)
3049       {
3050       vtkMedMesh* mesh=file->GetMesh(meshIndex);
3051       mesh->ClearMedSupports();
3052       }
3053
3054     int numProf = file->GetNumberOfProfile();
3055     for (int prof = 0; prof<numProf; prof++)
3056       {
3057       vtkMedProfile* profile = file->GetProfile(prof);
3058       if (profile->GetIds()!=NULL)
3059         profile->GetIds()->Initialize();
3060       }
3061     }
3062 }
3063
3064 void vtkMedReader::ClearMedFields()
3065 {
3066   this->Internal->FieldCache.clear();
3067   this->Internal->QuadOffsetKey.clear();
3068   this->Internal->QuadratureOffsetCache.clear();
3069
3070   std::map<std::string, vtkSmartPointer<vtkMedFile> >::iterator
3071       fieldfileit = this->Internal->MedFiles.begin();
3072   while(fieldfileit != this->Internal->MedFiles.end())
3073     {
3074     vtkMedFile* file = fieldfileit->second;
3075     fieldfileit++;
3076
3077     int numFields=file->GetNumberOfField();
3078     for(int ff=0; ff<numFields; ff++)
3079       {
3080       vtkMedField* field=file->GetField(ff);
3081       int nstep=field->GetNumberOfFieldStep();
3082       for(int sid=0; sid<nstep; sid++)
3083         {
3084         vtkMedFieldStep* step = field->GetFieldStep(sid);
3085         for(int id=0; id<step->GetNumberOfFieldOverEntity(); id++)
3086           {
3087           vtkMedFieldOverEntity * fieldOverEntity=step->GetFieldOverEntity(id);
3088           for(int pid = 0; pid < fieldOverEntity->GetNumberOfFieldOnProfile(); pid++)
3089             {
3090             vtkMedFieldOnProfile* fop = fieldOverEntity->GetFieldOnProfile(pid);
3091             if(fop->GetData() != NULL)
3092               fop->SetData(NULL);
3093             }
3094           }
3095         }
3096       }
3097     }
3098 }
3099
3100 void vtkMedReader::ClearCaches(int when)
3101 {
3102   switch(when)
3103   {
3104     case Initialize:
3105       this->Internal->CurrentDataSet.clear();
3106       this->Internal->DataSetCache.clear();
3107       this->Internal->FieldCache.clear();
3108       this->Internal->UsedSupports.clear();
3109       this->Internal->QuadratureOffsetCache.clear();
3110       this->Internal->QuadOffsetKey.clear();
3111       //this->Internal->Med2VTKPointIndex.clear();
3112       break;
3113     case StartRequest:
3114       this->Internal->CurrentDataSet.clear();
3115       this->Internal->UsedSupports.clear();
3116       if(this->CacheStrategy==CacheNothing)
3117         {
3118         this->ClearMedSupports();
3119         this->ClearMedFields();
3120         }
3121       else if(this->CacheStrategy==CacheGeometry)
3122         {
3123         this->ClearMedFields();
3124         }
3125       break;
3126     case AfterCreateMedSupports:
3127       // TODO : clear the unused supports and associated cached datasets and fields
3128       break;
3129     case EndBuildVTKSupports:
3130       break;
3131     case EndRequest:
3132       if(this->CacheStrategy==CacheNothing)
3133         {
3134         this->ClearMedSupports();
3135         this->ClearMedFields();
3136         }
3137       else if(this->CacheStrategy==CacheGeometry && this->AnimationMode != Modes)
3138         {
3139         this->ClearMedFields();
3140         }
3141       break;
3142   }
3143 }
3144
3145 void vtkMedReader::PrintSelf(ostream& os, vtkIndent indent)
3146 {
3147   PRINT_STRING(os, indent, FileName);
3148   PRINT_IVAR(os, indent, AnimationMode);
3149   PRINT_IVAR(os, indent, TimeIndexForIterations);
3150   PRINT_OBJECT(os, indent, PointFields);
3151   PRINT_OBJECT(os, indent, CellFields);
3152   PRINT_OBJECT(os, indent, QuadratureFields);
3153   PRINT_OBJECT(os, indent, ElnoFields);
3154   PRINT_OBJECT(os, indent, Groups);
3155   PRINT_OBJECT(os, indent, Entities);
3156   PRINT_IVAR(os, indent, CacheStrategy);
3157   this->Superclass::PrintSelf(os, indent);
3158 }