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