src/Plugins/MedReader/IO/vtkMedReader.cxx

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