src/Plugins/MedReader/IO/vtkMedDriver30.cxx

   1 // Copyright (C) 2010-2011  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 "vtkMedDriver30.h"
  21
  22 #include "vtkObjectFactory.h"
  23 #include "vtkStringArray.h"
  24 #include "vtkDataArray.h"
  25 #include "vtkIdTypeArray.h"
  26 #include "vtkDoubleArray.h"
  27 #include "vtkIntArray.h"
  28 #include "vtkCharArray.h"
  29
  30 #include "vtkMedFile.h"
  31 #include "vtkMedCartesianGrid.h"
  32 #include "vtkMedPolarGrid.h"
  33 #include "vtkMedCurvilinearGrid.h"
  34 #include "vtkMedUnstructuredGrid.h"
  35 #include "vtkMedField.h"
  36 #include "vtkMedMesh.h"
  37 #include "vtkMedFamily.h"
  38 #include "vtkMedUtilities.h"
  39 #include "vtkMedEntityArray.h"
  40 #include "vtkMedLocalization.h"
  41 #include "vtkMedProfile.h"
  42 #include "vtkMedFieldOverEntity.h"
  43 #include "vtkMedFieldStep.h"
  44 #include "vtkMedGroup.h"
  45 #include "vtkMedIntArray.h"
  46 #include "vtkMedInterpolation.h"
  47 #include "vtkMedFieldOnProfile.h"
  48 #include "vtkMedFraction.h"
  49 #include "vtkMedLink.h"
  50 #include "vtkMedStructElement.h"
  51 #include "vtkMedConstantAttribute.h"
  52 #include "vtkMedVariableAttribute.h"
  53
  54 #include "vtkMultiProcessController.h"
  55
  56 #include <string>
  57 using namespace std;
  58
  59 // vtkCxxRevisionMacro(vtkMedDriver30, "$Revision$")
  60 vtkStandardNewMacro(vtkMedDriver30)
  61
  62 vtkMedDriver30::vtkMedDriver30() : vtkMedDriver()
  63 {
  64 }
  65
  66 vtkMedDriver30::~vtkMedDriver30()
  67 {
  68 }
  69
  70 // Description:
  71 // load all Information data associated with this cartesian grid.
  72 void vtkMedDriver30::ReadRegularGridInformation(vtkMedRegularGrid* grid)
  73 {
  74   FileOpen open(this);
  75
  76   grid->SetDimension(grid->GetParentMesh()->GetNumberOfAxis());
  77
  78   for(int axis=0; axis < grid->GetDimension(); axis++)
  79     {
  80     med_int size;
  81     med_bool coordinatechangement, geotransformation;
  82
  83     if ((size = MEDmeshnEntity(
  84                   this->FileId,
  85                   grid->GetParentMesh()->GetName(),
  86                   grid->GetComputeStep().TimeIt,
  87                   grid->GetComputeStep().IterationIt,
  88                   MED_NODE,
  89                   MED_NONE,
  90                   (med_data_type)(MED_COORDINATE_AXIS1 + axis),
  91                   MED_NO_CMODE,
  92                   &coordinatechangement,
  93                   &geotransformation)) < 0)
  94       {
  95       vtkErrorMacro("ERROR : number of coordinates on X axis ...");
  96       }
  97
  98     grid->SetAxisSize(axis, size);
  99     }
 100
 101   med_int ncell = 1;
 102   if(grid->GetAxisSize(0) > 1)
 103     ncell *= grid->GetAxisSize(0)-1;
 104   if(grid->GetAxisSize(1) > 1)
 105     ncell *= grid->GetAxisSize(1)-1;
 106   if(grid->GetAxisSize(2) > 1)
 107     ncell *= grid->GetAxisSize(2)-1;
 108
 109   vtkMedEntity entity;
 110   entity.EntityType = MED_CELL;
 111
 112   switch(grid->GetDimension())
 113     {
 114     case 0 :
 115       entity.GeometryType = MED_POINT1;
 116       break;
 117     case 1 :
 118       entity.GeometryType = MED_SEG2;
 119       break;
 120     case 2 :
 121       entity.GeometryType = MED_QUAD4;
 122       break;
 123     case 3 :
 124       entity.GeometryType = MED_HEXA8;
 125       break;
 126     default :
 127         vtkErrorMacro("Unsupported dimension for curvilinear grid : "
 128                       << grid->GetDimension());
 129     return;
 130     }
 131
 132   vtkMedEntityArray* array = vtkMedEntityArray::New();
 133   array->SetParentGrid(grid);
 134   array->SetNumberOfEntity(ncell);
 135   array->SetEntity(entity);
 136   array->SetConnectivity(MED_NODAL);
 137   grid->AppendEntityArray(array);
 138   array->Delete();
 139   // this triggers the creation of undefined families
 140   this->LoadFamilyIds(array);
 141
 142 }
 143
 144 void  vtkMedDriver30::LoadRegularGridCoordinates(vtkMedRegularGrid* grid)
 145 {
 146   FileOpen open(this);
 147
 148   for(int axis=0; axis < grid->GetParentMesh()->GetNumberOfAxis(); axis++)
 149     {
 150
 151     vtkDataArray* coords = vtkMedUtilities::NewCoordArray();
 152     grid->SetAxisCoordinate(axis, coords);
 153     coords->Delete();
 154
 155     coords->SetNumberOfComponents(1);
 156     coords->SetNumberOfTuples(grid->GetAxisSize(axis));
 157
 158     if (MEDmeshGridIndexCoordinateRd(
 159           this->FileId,
 160           grid->GetParentMesh()->GetName(),
 161           grid->GetComputeStep().TimeIt,
 162           grid->GetComputeStep().IterationIt,
 163           axis+1,
 164           (med_float*)coords->GetVoidPointer(0)) < 0)
 165       {
 166       vtkErrorMacro("ERROR : read axis " << axis << " coordinates ...");
 167       grid->SetAxisCoordinate(axis, NULL);
 168       return;
 169       }
 170     }
 171
 172 }
 173
 174 // Description:
 175 // load all Information data associated with this standard grid.
 176 void vtkMedDriver30::ReadCurvilinearGridInformation(vtkMedCurvilinearGrid* grid)
 177 {
 178   FileOpen open(this);
 179
 180   grid->SetDimension(grid->GetParentMesh()->GetNumberOfAxis());
 181
 182   med_int size;
 183   med_bool coordinatechangement, geotransformation;
 184
 185   if ((size = MEDmeshnEntity(
 186                 this->FileId,
 187                 grid->GetParentMesh()->GetName(),
 188                 grid->GetComputeStep().TimeIt,
 189                 grid->GetComputeStep().IterationIt,
 190                 MED_NODE,
 191                 MED_NONE,
 192                 MED_COORDINATE,
 193                 MED_NO_CMODE,
 194                 &coordinatechangement,
 195                 &geotransformation)) < 0)
 196     {
 197     vtkErrorMacro("ReadCurvilinearGridInformation MEDmeshnEntity");
 198     }
 199
 200   grid->SetNumberOfPoints(size);
 201
 202   med_int axissize[3];
 203
 204   if(MEDmeshGridStructRd(
 205       this->FileId,
 206       grid->GetParentMesh()->GetName(),
 207       grid->GetComputeStep().TimeIt,
 208       grid->GetComputeStep().IterationIt,
 209       axissize) < 0)
 210     {
 211     vtkErrorMacro("ReadCurvilinearGridInformation MEDmeshGridStructRd");
 212     }
 213
 214   switch(grid->GetDimension())
 215     {
 216     case 0 :
 217       break;
 218     case 1 :
 219       grid->SetAxisSize(0, (axissize[0] <= 0 ? 1: axissize[0]));
 220       break;
 221     case 2 :
 222       grid->SetAxisSize(0, (axissize[0] <= 0 ? 1: axissize[0]));
 223       grid->SetAxisSize(1, (axissize[1] <= 0 ? 1: axissize[1]));
 224       break;
 225     case 3 :
 226       grid->SetAxisSize(0, (axissize[0] <= 0 ? 1: axissize[0]));
 227       grid->SetAxisSize(1, (axissize[1] <= 0 ? 1: axissize[1]));
 228       grid->SetAxisSize(2, (axissize[2] <= 0 ? 1: axissize[2]));
 229       break;
 230     default :
 231         vtkErrorMacro("Unsupported dimension for curvilinear grid : "
 232                       << grid->GetDimension());
 233     return;
 234     }
 235
 236   // A test to verify the number of points : total number of
 237   // points should be equal to the product of each axis size
 238   med_int size2 = 0;
 239
 240   if (grid->GetDimension() == 1)
 241     {
 242     size2 = grid->GetAxisSize(0);
 243     }
 244
 245   if (grid->GetDimension() == 2)
 246     {
 247     size2 = grid->GetAxisSize(0)*grid->GetAxisSize(1);
 248     }
 249
 250   if (grid->GetDimension() == 3)
 251     {
 252     size2 = grid->GetAxisSize(0)*grid->GetAxisSize(1)*grid->GetAxisSize(2);
 253     }
 254
 255   if(size != size2)
 256     {
 257     vtkErrorMacro("The total number of points of a Curvilinear grid should "
 258                   << "be the product of each axis size!");
 259     }
 260
 261   med_int ncell = 1;
 262   if ((grid->GetDimension() >= 1) && (grid->GetAxisSize(0) > 1))
 263     ncell *= grid->GetAxisSize(0)-1;
 264   if ((grid->GetDimension() >= 2) && (grid->GetAxisSize(1) > 1))
 265     ncell *= grid->GetAxisSize(1)-1;
 266   if ((grid->GetDimension() >= 3) && (grid->GetAxisSize(2) > 1))
 267     ncell *= grid->GetAxisSize(2)-1;
 268
 269   vtkMedEntity entity;
 270   entity.EntityType = MED_CELL;
 271
 272   switch(grid->GetDimension())
 273     {
 274     case 0 :
 275       entity.GeometryType = MED_POINT1;
 276       break;
 277     case 1 :
 278       entity.GeometryType = MED_SEG2;
 279       break;
 280     case 2 :
 281       entity.GeometryType = MED_QUAD4;
 282       break;
 283     case 3 :
 284       entity.GeometryType = MED_HEXA8;
 285       break;
 286     default :
 287         vtkErrorMacro("Unsupported dimension for curvilinear grid : "
 288                       << grid->GetDimension());
 289     return;
 290     }
 291
 292   vtkMedEntityArray* array = vtkMedEntityArray::New();
 293   array->SetParentGrid(grid);
 294   array->SetNumberOfEntity(ncell);
 295   array->SetEntity(entity);
 296   array->SetConnectivity(MED_NODAL);
 297   grid->AppendEntityArray(array);
 298   array->Delete();
 299   // this triggers the creation of undefined families
 300   this->LoadFamilyIds(array);
 301 }
 302
 303 // Description : read the number of entity of all geometry type
 304 // for a given entity type and a given connectivity mode
 305 void vtkMedDriver30::ReadNumberOfEntity(
 306     vtkMedUnstructuredGrid* grid,
 307     med_entity_type entityType,
 308     med_connectivity_mode connectivity)
 309 {
 310   FileOpen open(this);
 311
 312   med_bool changement, transformation;
 313
 314   const char* meshName = grid->GetParentMesh()->GetName();
 315
 316   const vtkMedComputeStep& cs = grid->GetComputeStep();
 317
 318   med_int nentity = MEDmeshnEntity(
 319                         this->FileId,
 320                         meshName,
 321                         cs.TimeIt,
 322                         cs.IterationIt,
 323                         entityType,
 324                         MED_GEO_ALL,
 325                         MED_UNDEF_DATATYPE ,
 326                         connectivity,
 327                         &changement,
 328                         &transformation );
 329
 330   for(med_int geotypeit = 1; geotypeit <= nentity; geotypeit++)
 331     {
 332     // read cell informations
 333     vtkMedEntity entity;
 334     entity.EntityType = entityType;
 335
 336     char geometryName[MED_NAME_SIZE+1] = "";
 337
 338     // this gives us the med_geometry_type
 339     if( MEDmeshEntityInfo( FileId, meshName,
 340                            cs.TimeIt,
 341                            cs.IterationIt,
 342                            entityType,
 343                            geotypeit,
 344                            geometryName,
 345                            &entity.GeometryType) < 0)
 346       {
 347       vtkErrorMacro("MEDmeshEntityInfo");
 348       continue;
 349       }
 350
 351     entity.GeometryName = geometryName;
 352     med_int ncell = 0;
 353
 354     if(entity.GeometryType == MED_POLYGON)
 355       {
 356       // read the number of cells of this type
 357       ncell = MEDmeshnEntity(this->FileId,
 358                              meshName,
 359                              cs.TimeIt,
 360                              cs.IterationIt,
 361                              entity.EntityType,
 362                              entity.GeometryType,
 363                              MED_INDEX_NODE,
 364                              connectivity,
 365                              &changement,
 366                              &transformation ) - 1;
 367       }
 368     else if(entity.GeometryType == MED_POLYHEDRON)
 369       {
 370       // read the number of cells of this type
 371       ncell = MEDmeshnEntity(this->FileId,
 372                              meshName,
 373                              cs.TimeIt,
 374                              cs.IterationIt,
 375                              entity.EntityType,
 376                              entity.GeometryType,
 377                              MED_INDEX_FACE,
 378                              connectivity,
 379                              &changement,
 380                              &transformation  ) - 1;
 381       }
 382     else
 383       {
 384       ncell = MEDmeshnEntity(this->FileId,
 385                              meshName,
 386                              cs.TimeIt,
 387                              cs.IterationIt,
 388                              entity.EntityType,
 389                              entity.GeometryType,
 390                              MED_CONNECTIVITY,
 391                              connectivity,
 392                              &changement,
 393                              &transformation  );
 394       }
 395
 396     if(ncell > 0)
 397       {
 398       vtkMedEntityArray* array = vtkMedEntityArray::New();
 399       array->SetParentGrid(grid);
 400       array->SetNumberOfEntity(ncell);
 401       array->SetEntity(entity);
 402       array->SetConnectivity(connectivity);
 403       grid->AppendEntityArray(array);
 404       array->Delete();
 405       // this triggers the creation of undefined families
 406       this->LoadFamilyIds(array);
 407       }
 408     }
 409 }
 410
 411 // Description:
 412 // load all Information data associated with this unstructured grid.
 413 void vtkMedDriver30::ReadUnstructuredGridInformation(
 414     vtkMedUnstructuredGrid* grid)
 415 {
 416   FileOpen open(this);
 417
 418   vtkMedMesh *mesh = grid->GetParentMesh();
 419
 420   const char *meshName = mesh->GetName();
 421   med_connectivity_mode connectivity;
 422
 423   med_bool changement;
 424   med_bool transformation;
 425   med_int profilesize;
 426
 427   char profilename[MED_NAME_SIZE+1];
 428   memset(profilename, '\0', MED_NAME_SIZE+1);
 429
 430   vtkMedComputeStep cs = grid->GetComputeStep();
 431
 432   // first check if we have points
 433   vtkIdType npoints = MEDmeshnEntityWithProfile(
 434                         this->FileId,
 435                         meshName,
 436                         cs.TimeIt,
 437                         cs.IterationIt,
 438                         MED_NODE,
 439                         MED_NONE,
 440                         MED_COORDINATE,
 441                         MED_NODAL,
 442                         MED_COMPACT_PFLMODE,
 443                         profilename,
 444                         &profilesize,
 445                         &changement,
 446                         &transformation);
 447
 448   if(npoints > 0)
 449     {
 450     grid->SetNumberOfPoints(npoints);
 451     }
 452   else
 453     {
 454     if(grid->GetPreviousGrid() == NULL)
 455       {
 456       vtkErrorMacro("No point and no previous grid");
 457       }
 458     grid->SetUsePreviousCoordinates(true);
 459     grid->SetNumberOfPoints(grid->GetPreviousGrid()->GetNumberOfPoints());
 460     return;
 461     }
 462
 463   this->ReadNumberOfEntity(grid, MED_CELL, MED_NODAL);
 464   this->ReadNumberOfEntity(grid, MED_CELL, MED_DESCENDING);
 465   this->ReadNumberOfEntity(grid, MED_DESCENDING_FACE, MED_NODAL);
 466   this->ReadNumberOfEntity(grid, MED_DESCENDING_FACE, MED_DESCENDING);
 467   this->ReadNumberOfEntity(grid, MED_DESCENDING_EDGE, MED_NODAL);
 468   this->ReadNumberOfEntity(grid, MED_DESCENDING_EDGE, MED_DESCENDING);
 469   this->ReadNumberOfEntity(grid, MED_STRUCT_ELEMENT, MED_NODAL);
 470   this->ReadNumberOfEntity(grid, MED_STRUCT_ELEMENT, MED_DESCENDING);
 471
 472   // create the point vtkMedEntityArray support
 473   vtkMedEntity entity;
 474   entity.EntityType = MED_NODE;
 475   entity.GeometryType = MED_POINT1;
 476   vtkMedEntityArray* pea = vtkMedEntityArray::New();
 477   pea->SetEntity(entity);
 478   pea->SetParentGrid(grid);
 479   pea->SetNumberOfEntity(grid->GetNumberOfPoints());
 480   grid->AppendEntityArray(pea);
 481   pea->Delete();
 482
 483   this->LoadFamilyIds(pea);
 484 }
 485
 486 void vtkMedDriver30::ReadFamilyInformation(vtkMedMesh* mesh, vtkMedFamily* family)
 487 {
 488   FileOpen open(this);
 489
 490   med_int familyid;
 491   med_int ngroup;
 492   char* groupNames = NULL;
 493   const  char* meshName = mesh->GetName();
 494
 495   ngroup = MEDnFamilyGroup(FileId, meshName, family->GetMedIterator());
 496
 497   bool has_no_group = false;
 498   if(ngroup <= 0)
 499     {
 500     if(ngroup < 0)
 501       {
 502       vtkErrorMacro("Error while reading the number of groups");
 503       }
 504     ngroup = 1;
 505     has_no_group = true;
 506     }
 507
 508   groupNames = new char[ngroup * MED_LNAME_SIZE + 1];
 509   memset(groupNames, '\0', ngroup * MED_LNAME_SIZE + 1);
 510
 511   // special case for files written by med < 3,
 512   // I have to use the 23 interface
 513   if(mesh->GetParentFile()->GetVersionMajor() < 3)
 514     {
 515     med_int *attributenumber = NULL;
 516     med_int *attributevalue = NULL;
 517     char *attributedes = NULL;
 518
 519     med_int nattr = MEDnFamily23Attribute(
 520                       this->FileId,
 521                       meshName,
 522                       family->GetMedIterator());
 523
 524     if(nattr < 0)
 525       {
 526       vtkErrorMacro("MEDnFamily23Attribute");
 527       }
 528
 529     if(nattr > 0)
 530       {
 531       attributenumber = new med_int[nattr];
 532       attributevalue = new med_int[nattr];
 533       attributedes = new char[nattr*MED_COMMENT_SIZE+1];
 534       memset(attributedes, '\0', nattr*MED_COMMENT_SIZE+1);
 535       }
 536
 537     char familyName[MED_NAME_SIZE+1] = "";
 538
 539     if(MEDfamily23Info (this->FileId,
 540                         meshName,
 541                         family->GetMedIterator(),
 542                         familyName,
 543                         attributenumber,
 544                         attributevalue,
 545                         attributedes,
 546                         &familyid,
 547                         groupNames ) < 0)
 548       {
 549       vtkDebugMacro("MEDfamily23Info");
 550       }
 551
 552     family->SetName(familyName);
 553
 554     if(attributenumber != NULL)
 555       delete[] attributenumber;
 556     if(attributevalue != NULL)
 557       delete[] attributevalue;
 558     if(attributedes != NULL)
 559       delete[] attributedes;
 560     }
 561   else
 562     {
 563     char familyName[MED_NAME_SIZE+1] = "";
 564     if(MEDfamilyInfo( this->FileId,
 565                       meshName,
 566                       family->GetMedIterator(),
 567                       familyName,
 568                       &familyid,
 569                       groupNames ) < 0)
 570       {
 571       vtkErrorMacro(
 572           "vtkMedDriver30::ReadInformation(vtkMedFamily* family)"
 573           << " cannot read family informations.");
 574       return;
 575       }
 576     family->SetName(familyName);
 577     }
 578
 579   family->SetId(familyid);
 580
 581   if(familyid <= 0)
 582     {
 583     family->SetPointOrCell(vtkMedUtilities::OnCell);
 584     mesh->AppendCellFamily(family);
 585     }
 586   else
 587     {
 588     family->SetPointOrCell(vtkMedUtilities::OnPoint);
 589     mesh->AppendPointFamily(family);
 590     }
 591
 592   family->AllocateNumberOfGroup(ngroup);
 593   // if there where no group, set the name to the default value
 594   if(has_no_group)
 595     {
 596     memcpy(groupNames, vtkMedUtilities::NoGroupName,
 597            strlen(vtkMedUtilities::NoGroupName));
 598     }
 599
 600   for(int index = 0; index < ngroup; index++)
 601     {
 602     char realGroupName[MED_LNAME_SIZE + 1];
 603     memset(realGroupName, '\0', MED_LNAME_SIZE + 1);
 604     memcpy(realGroupName, groupNames + index * MED_LNAME_SIZE,
 605         MED_LNAME_SIZE * sizeof(char));
 606     vtkMedGroup* group = mesh->GetOrCreateGroup(family->GetPointOrCell(),
 607         realGroupName);
 608
 609     family->SetGroup(index, group);
 610     }
 611
 612   delete[] groupNames;
 613
 614   if(familyid == 0)
 615     {
 616     vtkMedFamily* famzero = vtkMedFamily::New();
 617     mesh->AppendPointFamily(famzero);
 618     famzero->Delete();
 619
 620     famzero->SetName(family->GetName());
 621     famzero->SetMedIterator(family->GetMedIterator());
 622     famzero->SetId(family->GetId());
 623     famzero->SetPointOrCell(vtkMedUtilities::OnPoint);
 624     famzero->AllocateNumberOfGroup(family->GetNumberOfGroup());
 625     for(int gid=0; gid<family->GetNumberOfGroup(); gid++)
 626       {
 627       vtkMedGroup* group = mesh->GetOrCreateGroup(
 628           vtkMedUtilities::OnPoint,
 629           family->GetGroup(gid)->GetName());
 630       famzero->SetGroup(gid, group);
 631       mesh->AppendPointGroup(group);
 632       }
 633     }
 634 }
 635
 636 void  vtkMedDriver30::ReadLinkInformation(vtkMedLink* link)
 637 {
 638   med_int size;
 639   char linkMeshName[MED_NAME_SIZE+1] = "";
 640   if(MEDlinkInfo(this->FileId,
 641                  link->GetMedIterator(),
 642                  linkMeshName,
 643                  &size) < 0)
 644     {
 645     vtkErrorMacro("MEDlinkInfo");
 646     return;
 647     }
 648   link->SetMeshName(linkMeshName);
 649   if(size <= 0)
 650     return;
 651
 652   char* path = new char[size + 1];
 653   memset(path, '\0', size+1);
 654   if(MEDlinkRd(this->FileId, link->GetMeshName(), path) < 0)
 655     {
 656     vtkErrorMacro("MEDlinkRd");
 657     memset(path, '\0', size+1);
 658     }
 659
 660   link->SetLink(path);
 661
 662   delete[] path;
 663 }
 664
 665 void vtkMedDriver30::ReadFileInformation(vtkMedFile* file)
 666 {
 667   FileOpen open(this);
 668
 669   char comment[MED_COMMENT_SIZE+1] = "";
 670
 671   MEDfileCommentRd(this->FileId,
 672                   comment);
 673
 674   file->SetComment(comment);
 675
 676   med_int major, minor, release;
 677   MEDfileNumVersionRd(this->FileId, &major, &minor, &release);
 678   file->SetVersionMajor(major);
 679   file->SetVersionMinor(minor);
 680   file->SetVersionRelease(release);
 681
 682   int nlink = MEDnLink(this->FileId);
 683   file->AllocateNumberOfLink(nlink);
 684   for(int linkid=0; linkid<nlink; linkid++)
 685     {
 686     vtkMedLink* link = file->GetLink(linkid);
 687     link->SetMedIterator(linkid+1);
 688     this->ReadLinkInformation(link);
 689     }
 690
 691   int nprof = MEDnProfile(FileId);
 692   // Reading id s not possible in parallel if the file contains Profiles
 693   vtkMultiProcessController* controller = vtkMultiProcessController::GetGlobalController();
 694   if (controller != NULL)
 695     if ((nprof != 0) && (controller->GetNumberOfProcesses() > 1))
 696     {
 697       vtkWarningMacro("ATTENTION: The MED Reader cannot read profiles when used in parallel");
 698     return;
 699     }
 700   file->AllocateNumberOfProfile(nprof);
 701   for(int i = 0; i < nprof; i++)
 702     {
 703     vtkMedProfile* profile = file->GetProfile(i);
 704     profile->SetMedIterator(i + 1);
 705     profile->SetParentFile(file);
 706     this->ReadProfileInformation(profile);
 707     }
 708
 709   int nloc = MEDnLocalization(this->FileId);
 710   file->AllocateNumberOfLocalization(nloc);
 711   for(int i = 0; i < nloc; i++)
 712     {
 713     vtkMedLocalization* loc = file->GetLocalization(i);
 714     loc->SetMedIterator(i + 1);
 715     loc->SetParentFile(file);
 716     this->ReadLocalizationInformation(loc);
 717     }
 718
 719   int nsupportmesh = MEDnSupportMesh(this->FileId);
 720   file->AllocateNumberOfSupportMesh(nsupportmesh);
 721   for(int i = 0; i < nsupportmesh; i++)
 722     {
 723     vtkMedMesh* supportmesh = file->GetSupportMesh(i);
 724     supportmesh->SetMedIterator(i + 1);
 725     supportmesh->SetParentFile(file);
 726     this->ReadSupportMeshInformation(supportmesh);
 727     }
 728
 729   int nmesh = MEDnMesh(this->FileId);
 730   file->AllocateNumberOfMesh(nmesh);
 731   for(int i = 0; i < nmesh; i++)
 732     {
 733     vtkMedMesh* mesh = file->GetMesh(i);
 734     mesh->SetMedIterator(i + 1);
 735     mesh->SetParentFile(file);
 736     this->ReadMeshInformation(mesh);
 737     }
 738
 739   int nfields = MEDnField(this->FileId);
 740   file->AllocateNumberOfField(nfields);
 741   for(int i = 0; i < nfields; i++)
 742     {
 743     vtkMedField* field = file->GetField(i);
 744     field->SetMedIterator(i + 1);
 745     field->SetParentFile(file);
 746     this->ReadFieldInformation(field);
 747     field->ComputeFieldType();
 748     while(field->HasManyFieldTypes())
 749       {
 750       vtkMedField* newfield = vtkMedField::New();
 751       int type = field->GetFirstType();
 752       newfield->ExtractFieldType(field, type);
 753       file->AppendField(newfield);
 754       newfield->Delete();
 755       }
 756     }
 757
 758   int nstruct = MEDnStructElement(this->FileId);
 759
 760   file->AllocateNumberOfStructElement(nstruct);
 761   for(int i = 0; i < nstruct; i++)
 762     {
 763     vtkMedStructElement* structelem = file->GetStructElement(i);
 764     structelem->SetMedIterator(i+1);
 765     structelem->SetParentFile(file);
 766     this->ReadStructElementInformation(structelem);
 767     }
 768
 769 }
 770
 771 void  vtkMedDriver30::ReadStructElementInformation(
 772     vtkMedStructElement* structelem)
 773 {
 774
 775   FileOpen open(this);
 776
 777   char modelname[MED_NAME_SIZE+1] = "";
 778   med_geometry_type mgeotype;
 779   med_int modeldim;
 780   char supportmeshname[MED_NAME_SIZE+1] = "";
 781   med_entity_type sentitytype;
 782   med_int snbofnode;
 783   med_int snbofcell;
 784   med_geometry_type sgeotype;
 785   med_int nbofconstantattribute;
 786   med_bool anyprofile;
 787   med_int nbofvariableattribute;
 788
 789   if(MEDstructElementInfo (this->FileId,
 790                            structelem->GetMedIterator(),
 791                            modelname,
 792                            &mgeotype,
 793                            &modeldim,
 794                            supportmeshname,
 795                            &sentitytype,
 796                            &snbofnode,
 797                            &snbofcell,
 798                            &sgeotype,
 799                            &nbofconstantattribute,
 800                            &anyprofile,
 801                            &nbofvariableattribute) < 0)
 802     {
 803     vtkErrorMacro("Error in MEDstructElementInfo");
 804     return;
 805     }
 806   structelem->SetName(modelname);
 807   structelem->SetGeometryType(mgeotype);
 808   structelem->SetModelDimension(modeldim);
 809   structelem->SetSupportMeshName(supportmeshname);
 810   structelem->SetSupportEntityType(sentitytype);
 811   structelem->SetSupportNumberOfNode(snbofnode);
 812   structelem->SetSupportNumberOfCell(snbofcell);
 813   structelem->SetSupportGeometryType(sgeotype);
 814   structelem->AllocateNumberOfConstantAttribute(nbofconstantattribute);
 815   structelem->AllocateNumberOfVariableAttribute(nbofvariableattribute);
 816   structelem->SetAnyProfile(anyprofile);
 817
 818   for(int attit = 0; attit < nbofconstantattribute; attit ++)
 819     {
 820     vtkMedConstantAttribute* constatt = structelem->GetConstantAttribute(attit);
 821     constatt->SetMedIterator(attit+1);
 822     constatt->SetParentStructElement(structelem);
 823     this->ReadConstantAttributeInformation(constatt);
 824     }
 825
 826   for(int attit = 0; attit < nbofvariableattribute; attit ++)
 827     {
 828     vtkMedVariableAttribute* varatt = structelem->GetVariableAttribute(attit);
 829     varatt->SetMedIterator(attit+1);
 830     varatt->SetParentStructElement(structelem);
 831     this->ReadVariableAttributeInformation(varatt);
 832     }
 833 }
 834
 835 void vtkMedDriver30::ReadConstantAttributeInformation(
 836     vtkMedConstantAttribute* constAttr)
 837 {
 838
 839   FileOpen open(this);
 840
 841   char constattname[MED_NAME_SIZE+1] = "";
 842   med_attribute_type constatttype;
 843   med_int nbofcomponent;
 844   med_entity_type sentitytype;
 845   char profilename[MED_NAME_SIZE+1] = "";
 846   med_int profilesize;
 847
 848   if(MEDstructElementConstAttInfo(
 849       this->FileId,
 850       constAttr->GetParentStructElement()->GetName(),
 851       constAttr->GetMedIterator(),
 852       constattname,
 853       &constatttype,
 854       &nbofcomponent,
 855       &sentitytype,
 856       profilename,
 857       &profilesize)   < 0)
 858     {
 859     vtkErrorMacro("MEDstructElementConstAttInfo error");
 860     return;
 861     }
 862
 863   constAttr->SetName(constattname);
 864   constAttr->SetAttributeType(constatttype);
 865   constAttr->SetNumberOfComponent(nbofcomponent);
 866   constAttr->SetSupportEntityType(sentitytype);
 867   constAttr->SetProfileName(profilename);
 868   constAttr->SetProfileSize(profilesize);
 869
 870   vtkAbstractArray* values = vtkMedUtilities::NewArray(constatttype);
 871   if(values == NULL)
 872     return;
 873   constAttr->SetValues(values);
 874   values->Delete();
 875
 876   values->SetNumberOfComponents(nbofcomponent);
 877   vtkIdType ntuple = 0;
 878   if((strcmp(profilename, MED_NO_PROFILE) != 0) &&
 879      (strcmp(profilename, "\0") != 0))
 880     {
 881     ntuple = profilesize;
 882     }
 883   else if(constAttr->GetSupportEntityType() == MED_CELL)
 884     {
 885     ntuple = constAttr->GetParentStructElement()->GetSupportNumberOfCell();
 886     }
 887   else
 888     {
 889     ntuple = constAttr->GetParentStructElement()->GetSupportNumberOfNode();
 890     }
 891   values->SetNumberOfTuples(ntuple);
 892
 893   void* ptr = NULL;
 894   vtkSmartPointer<vtkCharArray> buffer = vtkSmartPointer<vtkCharArray>::New();
 895   if(constatttype != MED_ATT_NAME)
 896     {
 897     ptr = values->GetVoidPointer(0);
 898     }
 899   else
 900     {
 901     buffer->SetNumberOfValues(MED_NAME_SIZE*nbofcomponent*ntuple);
 902     ptr = buffer->GetVoidPointer(0);
 903     }
 904
 905   if(MEDstructElementConstAttRd (this->FileId,
 906         constAttr->GetParentStructElement()->GetName(),
 907         constAttr->GetName(), ptr) < 0)
 908     {
 909     vtkErrorMacro("MEDstructElementConstAttRd");
 910     return;
 911     }
 912
 913   if(constatttype == MED_ATT_NAME)
 914     {
 915     char name[MED_NAME_SIZE+1] = "";
 916     char* nameptr = (char*) ptr;
 917     vtkStringArray* names = vtkStringArray::SafeDownCast(values);
 918     for(vtkIdType id = 0; id < nbofcomponent*ntuple; id++)
 919       {
 920       memset(name, '\0', MED_NAME_SIZE+1);
 921       strncpy(name, nameptr + id * MED_NAME_SIZE, MED_NAME_SIZE);
 922       names->SetValue(id, name);
 923       }
 924     }
 925
 926   return;
 927 }
 928
 929 void vtkMedDriver30::ReadVariableAttributeInformation(
 930     vtkMedVariableAttribute* varAttr)
 931 {
 932
 933   FileOpen open(this);
 934
 935   char varattname[MED_NAME_SIZE+1] = "";
 936   med_attribute_type varatttype;
 937   med_int nbofcomponent;
 938
 939   if(MEDstructElementVarAttInfo (
 940       this->FileId,
 941       varAttr->GetParentStructElement()->GetName(),
 942       varAttr->GetMedIterator(),
 943       varattname,
 944       &varatttype,
 945       &nbofcomponent) < 0)
 946     {
 947     vtkErrorMacro("MEDstructElementVarAttInfo");
 948     return;
 949     }
 950
 951   varAttr->SetName(varattname);
 952   varAttr->SetAttributeType(varatttype);
 953   varAttr->SetNumberOfComponent(nbofcomponent);
 954
 955   return;
 956 }
 957
 958 void vtkMedDriver30::ReadProfileInformation(vtkMedProfile* profile)
 959 {
 960   FileOpen open(this);
 961
 962   med_int nelem;
 963   char profileName[MED_NAME_SIZE+1] = "";
 964
 965   if(MEDprofileInfo(this->FileId,
 966                     profile->GetMedIterator(),
 967                     profileName,
 968                     &nelem) < 0)
 969     {
 970     vtkErrorMacro("cannot read information on profile"
 971         << profile->GetMedIterator());
 972     }
 973   profile->SetName(profileName);
 974   profile->SetNumberOfElement(nelem);
 975 }
 976
 977 void vtkMedDriver30::ReadFieldInformation(vtkMedField* field)
 978 {
 979   FileOpen open(this);
 980
 981   if (field->GetMedIterator() == 0)
 982     return;
 983
 984   int ncomp = MEDfieldnComponent(FileId, field->GetMedIterator());
 985
 986   if(ncomp < 0)
 987     {
 988     field->SetNumberOfComponent(-1);
 989     vtkErrorMacro("cannot read the number of component of field "
 990         << field->GetMedIterator())
 991     return;
 992     }
 993
 994   field->SetNumberOfComponent(ncomp);
 995
 996   char* units = new char[MED_SNAME_SIZE * ncomp + 1];
 997   char* componentNames = new char[MED_SNAME_SIZE * ncomp + 1];
 998   memset(units, '\0', MED_SNAME_SIZE * ncomp + 1);
 999   memset(componentNames, '\0', MED_SNAME_SIZE * ncomp + 1);
1000
1001   //med_type_champ dataType;
1002   med_field_type dataType;
1003   med_int nstep;
1004   med_bool localmesh;
1005
1006   char name[MED_NAME_SIZE+1] = "";
1007   char meshName[MED_NAME_SIZE+1] = "";
1008   char unit[MED_SNAME_SIZE+1] = "";
1009
1010   if( MEDfieldInfo( FileId,
1011                     field->GetMedIterator(),
1012                     name,
1013                     meshName,
1014                     &localmesh,
1015                     &dataType,
1016                     componentNames,
1017                     units,
1018                     unit,
1019                     &nstep) < 0)
1020     {
1021     vtkErrorMacro("cannot read the informations on field "
1022         << field->GetMedIterator())
1023     return;
1024     }
1025
1026   field->SetName(name);
1027   field->SetMeshName(meshName);
1028   field->SetTimeUnit(unit);
1029   field->SetDataType(dataType);
1030   field->SetLocal(localmesh);
1031
1032   for(int comp = 0; comp < ncomp; comp++)
1033     {
1034     char unit[MED_NAME_SIZE + 1] = "";
1035     memcpy(unit, units + MED_SNAME_SIZE * comp, MED_SNAME_SIZE * sizeof(char));
1036     field->GetUnit()->SetValue(comp, unit);
1037
1038     char compName[MED_SNAME_SIZE + 1] = "";
1039     memcpy(compName, componentNames + MED_SNAME_SIZE * comp, MED_SNAME_SIZE
1040         * sizeof(char));
1041     field->GetComponentName()->SetValue(comp, compName);
1042     }
1043
1044   delete[] units;
1045   delete[] componentNames;
1046
1047   med_int ninterp = MEDfieldnInterp(FileId, field->GetName());
1048   if(ninterp < 0)
1049     {
1050     vtkErrorMacro("Error in MEDfieldnInterp");
1051     return;
1052     }
1053
1054   field->AllocateNumberOfInterpolation(ninterp);
1055
1056   for(med_int interpit=0; interpit<ninterp; interpit++)
1057     {
1058     vtkMedInterpolation* interp = field->GetInterpolation(interpit);
1059     interp->SetMedIterator(interpit + 1);
1060     this->ReadInterpolationInformation(interp);
1061     }
1062
1063   vtkMedFieldStep* previousStep = NULL;
1064
1065   for(med_int csit = 0; csit < nstep; csit++)
1066     {
1067     vtkMedFieldStep* step = vtkMedFieldStep::New();
1068     step->SetMedIterator(csit + 1);
1069     step->SetParentField(field);
1070     this->ReadFieldStepInformation(step, csit == 0);
1071     field->AddFieldStep(step);
1072     step->SetPreviousStep(previousStep);
1073     previousStep = step;
1074     step->Delete();
1075     }
1076 }
1077
1078 void vtkMedDriver30::ReadFieldStepInformation(vtkMedFieldStep* step, bool readAllEntityInfo)
1079 {
1080   vtkMedComputeStep cs;
1081   vtkMedComputeStep meshcs;
1082   vtkMedField* field = step->GetParentField();
1083
1084   FileOpen open(this);
1085
1086   if( MEDfieldComputingStepMeshInfo(
1087         FileId,
1088         field->GetName(),
1089         step->GetMedIterator(),
1090         &cs.TimeIt,
1091         &cs.IterationIt,
1092         &cs.TimeOrFrequency,
1093         &meshcs.TimeIt,
1094         &meshcs.IterationIt) < 0)
1095     {
1096     vtkErrorMacro("Error in MEDfieldComputingStepMeshInfo");
1097     return;
1098     }
1099
1100   step->SetComputeStep(cs);
1101   step->SetMeshComputeStep(meshcs);
1102
1103   if(!readAllEntityInfo || step->GetEntityInfoLoaded())
1104     return;
1105
1106   step->SetEntityInfoLoaded(1);
1107
1108   vtkMedFile* file = field->GetParentFile();
1109   vtkMedMesh* mesh = file->GetMesh(field->GetMeshName());
1110
1111   if(mesh == NULL)
1112     return;
1113
1114   std::set<vtkMedEntity> entities;
1115   mesh->GatherMedEntities(entities);
1116
1117   vtkMedEntity pointEntity;
1118   pointEntity.EntityType = MED_NODE;
1119   pointEntity.GeometryType = MED_NONE;
1120   pointEntity.GeometryName = MED_NAME_BLANK;
1121   entities.insert(pointEntity);
1122
1123   std::set<vtkMedEntity>::iterator entity_it = entities.begin();
1124   while(entity_it != entities.end())
1125     {
1126     vtkMedEntity entity = *entity_it;
1127     entity_it++;
1128
1129     med_int nvalues = 0;
1130     med_int nprofile;
1131     char profilename[MED_NAME_SIZE+1] = "";
1132     char localizationname[MED_NAME_SIZE+1] = "";
1133
1134     nprofile = MEDfieldnProfile(
1135         this->FileId,
1136         field->GetName(),
1137         step->GetComputeStep().TimeIt,
1138         step->GetComputeStep().IterationIt,
1139         entity.EntityType,
1140         entity.GeometryType,
1141         profilename,
1142         localizationname);
1143     if(nprofile < 0)
1144       {
1145       vtkErrorMacro("MEDfieldnProfile");
1146       continue;
1147       }
1148
1149     bool hasprofile = (nprofile > 0);
1150     if(!hasprofile)
1151       {
1152       nprofile = 1;
1153       }
1154
1155     med_int profilesize;
1156     med_int nintegrationpoint;
1157
1158     for(int pid=0; pid<nprofile; pid++)
1159       {
1160       med_int medid = (hasprofile ? pid+1 : -1);
1161       nvalues = MEDfieldnValueWithProfile(
1162                   this->FileId,
1163                   field->GetName(),
1164                   step->GetComputeStep().TimeIt,
1165                   step->GetComputeStep().IterationIt,
1166                   entity.EntityType,
1167                   entity.GeometryType,
1168                   medid,
1169                   MED_COMPACT_PFLMODE,
1170                   profilename,
1171                   &profilesize,
1172                   localizationname,
1173                   &nintegrationpoint);
1174
1175       if(nvalues < 0)
1176         {
1177         vtkErrorMacro("MEDfieldnValueWithProfile");
1178         continue;
1179         }
1180       else if(nvalues > 0)
1181         {
1182         // I have found a profile with values, stop the loop here
1183         break;
1184         }
1185       }
1186
1187     if(nvalues > 0)
1188       {
1189       vtkMedFieldOverEntity* fieldOverEntity = vtkMedFieldOverEntity::New();
1190       step->AppendFieldOverEntity(fieldOverEntity);
1191       fieldOverEntity->Delete();
1192
1193       fieldOverEntity->SetParentStep(step);
1194       fieldOverEntity->SetEntity(entity);
1195
1196       this->ReadFieldOverEntityInformation(fieldOverEntity);
1197       }
1198     }
1199 }
1200
1201 void vtkMedDriver30::ReadFieldOverEntityInformation(vtkMedFieldOverEntity* fieldOverEntity)
1202 {
1203   FileOpen open(this);
1204
1205   vtkMedFieldStep* step = fieldOverEntity->GetParentStep();
1206   vtkMedField* field = step->GetParentField();
1207   vtkMedEntity entity = fieldOverEntity->GetEntity();
1208
1209   const char* fieldName = field->GetName();
1210   const vtkMedComputeStep& cs = step->GetComputeStep();
1211
1212   char profilename[MED_NAME_SIZE+1] = "";
1213   char localizationname[MED_NAME_SIZE+1] = "";
1214
1215   med_int nProfiles = MEDfieldnProfile(
1216                         this->FileId,
1217                         fieldName,
1218                         cs.TimeIt,
1219                         cs.IterationIt,
1220                         entity.EntityType,
1221                         entity.GeometryType,
1222                         profilename,
1223                         localizationname);
1224
1225   if(nProfiles < 0)
1226     {
1227     vtkErrorMacro("MEDfieldnProfile");
1228     }
1229   else if(nProfiles == 0)
1230     {
1231     fieldOverEntity->SetHasProfile(0);
1232     nProfiles = 1;
1233     }
1234   else
1235     {
1236     fieldOverEntity->SetHasProfile(1);
1237     }
1238   fieldOverEntity->AllocateNumberOfFieldOnProfile(nProfiles);
1239   for(int profit = 0; profit < nProfiles; profit++)
1240     {
1241     vtkMedFieldOnProfile* fop = fieldOverEntity->GetFieldOnProfile(profit);
1242     med_int medid = (fieldOverEntity->GetHasProfile()? profit+1: -1);
1243     fop->SetMedIterator(medid);
1244     fop->SetParentFieldOverEntity(fieldOverEntity);
1245     this->ReadFieldOnProfileInformation(fop);
1246     }
1247 }
1248
1249 void vtkMedDriver30::ReadFieldOnProfileInformation(vtkMedFieldOnProfile* fop)
1250 {
1251   vtkMedFieldOverEntity* fieldOverEntity = fop->GetParentFieldOverEntity();
1252   vtkMedFieldStep* step = fieldOverEntity->GetParentStep();
1253   vtkMedField* field = step->GetParentField();
1254
1255   const vtkMedComputeStep& cs = step->GetComputeStep();
1256   med_int profilesize;
1257   med_int nbofintegrationpoint;
1258
1259   char profileName[MED_NAME_SIZE+1] = "";
1260   char localizationName[MED_NAME_SIZE+1] = "";
1261
1262   med_int nvalue = MEDfieldnValueWithProfile(FileId,
1263                     field->GetName(),
1264                     cs.TimeIt,
1265                     cs.IterationIt,
1266                     fieldOverEntity->GetEntity().EntityType,
1267                     fieldOverEntity->GetEntity().GeometryType,
1268                     fop->GetMedIterator(),
1269                     MED_COMPACT_STMODE,
1270                     profileName,
1271                     &profilesize,
1272                     localizationName,
1273                     &nbofintegrationpoint);
1274
1275   if(nvalue < 0)
1276     {
1277     vtkErrorMacro("Error while reading MEDfieldnValueWithProfile");
1278     }
1279
1280   fop->SetProfileName(profileName);
1281   fop->SetLocalizationName(localizationName);
1282   fop->SetNumberOfValues(nvalue);
1283   fop->SetNumberOfIntegrationPoint(nbofintegrationpoint);
1284   fop->SetProfileSize(profilesize);
1285 }
1286
1287 void vtkMedDriver30::ReadMeshInformation(vtkMedMesh* mesh)
1288 {
1289   FileOpen open(this);
1290
1291   med_int mdim = 0;
1292   med_int sdim = 0;
1293   med_mesh_type meshtype;
1294
1295   med_sorting_type sortingtype;
1296   med_int nstep = 0;
1297   med_axis_type axistype;
1298   med_int naxis;
1299
1300   if ( (naxis=MEDmeshnAxis(this->FileId, mesh->GetMedIterator())) <0 )
1301     {
1302     vtkDebugMacro("Error reading mesh axis number");
1303     }
1304
1305   if(naxis == 0)
1306     naxis=MEDmeshnAxis(this->FileId, mesh->GetMedIterator());
1307
1308   char axisname[3*MED_SNAME_SIZE+1] = "";
1309   char axisunit[3*MED_SNAME_SIZE+1] = "";
1310   char name[MED_NAME_SIZE+1] = "";
1311   char description[MED_COMMENT_SIZE+1] = "";
1312   char timeUnit[MED_SNAME_SIZE+1] = "";
1313
1314   if( MEDmeshInfo( this->FileId,
1315         mesh->GetMedIterator(),
1316         name,
1317         &sdim,
1318         &mdim,
1319         &meshtype,
1320         description,
1321         timeUnit,
1322         &sortingtype,
1323         &nstep,
1324         &axistype,
1325         axisname,
1326         axisunit ) )
1327     {
1328     vtkErrorMacro("Error reading mesh");
1329     }
1330
1331   mesh->SetName(name);
1332   mesh->SetDescription(description);
1333   mesh->SetTimeUnit(timeUnit);
1334   mesh->SetSpaceDimension(sdim);
1335   mesh->SetMeshDimension(mdim);
1336   mesh->SetMeshType(meshtype);
1337   mesh->SetSortingType(sortingtype);
1338   mesh->SetAxisType(axistype);
1339   mesh->SetNumberOfAxis(naxis);
1340
1341   for(int axis = 0; axis < naxis; axis++)
1342     {
1343     char theaxisname[MED_SNAME_SIZE+1] = "";
1344     char theaxisunit[MED_SNAME_SIZE+1] = "";
1345     strncpy(theaxisname, axisname + axis*MED_SNAME_SIZE, MED_SNAME_SIZE);
1346     strncpy(theaxisunit, axisunit + axis*MED_SNAME_SIZE, MED_SNAME_SIZE);
1347     mesh->GetAxisName()->SetValue(axis, theaxisname);
1348     mesh->GetAxisUnit()->SetValue(axis, theaxisunit);
1349     }
1350
1351   char universalName[MED_LNAME_SIZE+1] = "";
1352
1353   if(MEDmeshUniversalNameRd(this->FileId, name,
1354       universalName) < 0)
1355     {
1356     vtkDebugMacro("MEDmeshUniversalNameRd < 0");
1357     }
1358   mesh->SetUniversalName(universalName);
1359
1360   // read the Information data of all families.
1361   // writing the family 0 is optional,
1362   // but I need it, so add it if it is not here.
1363
1364   med_int nfam = MEDnFamily(this->FileId, name);
1365
1366   for(int index = 0; index < nfam; index++)
1367     {
1368     vtkMedFamily* family = vtkMedFamily::New();
1369     family->SetMedIterator(index + 1);
1370     this->ReadFamilyInformation(mesh, family);
1371     family->Delete();
1372     }
1373
1374   // this creates a family 0 if none has been read
1375   vtkMedFamily* familyZeroOnCell = mesh->GetOrCreateCellFamilyById(0);
1376   vtkMedFamily* familyZeroOnPoint = mesh->GetOrCreatePointFamilyById(0);
1377
1378   // Load Information regarding the grid type
1379   if(meshtype == MED_STRUCTURED_MESH)
1380     {
1381     // Do it for structured data
1382     med_grid_type mtg;
1383     if(MEDmeshGridTypeRd(FileId, name, &mtg) < 0)
1384       {
1385       vtkErrorMacro("Error during structured grid Information loading.");
1386       return;
1387       }
1388     mesh->SetStructuredGridType(mtg);
1389     }
1390
1391   vtkMedGrid* previousGrid = NULL;
1392   for(int gid=1; gid <= nstep; gid++)
1393     {
1394     vtkMedComputeStep cs;
1395     if(MEDmeshComputationStepInfo(FileId,
1396                                   name,
1397                                   gid,
1398                                   &cs.TimeIt,
1399                                   &cs.IterationIt,
1400                                   &cs.TimeOrFrequency) < 0)
1401       {
1402       vtkErrorMacro("MEDmeshComputationStepInfo error");
1403       }
1404     // Load Information regarding the grid type
1405     vtkMedGrid* grid = NULL;
1406     if(meshtype == MED_STRUCTURED_MESH)
1407       {
1408       switch(mesh->GetStructuredGridType())
1409         {
1410         case MED_CARTESIAN_GRID:
1411           grid = vtkMedCartesianGrid::New();
1412           break;
1413         case MED_POLAR_GRID:
1414           grid = vtkMedPolarGrid::New();
1415           break;
1416         case MED_CURVILINEAR_GRID:
1417           grid = vtkMedCurvilinearGrid::New();
1418           break;
1419         default:
1420           vtkErrorMacro("Unknown structured grid type " << mesh->GetStructuredGridType());
1421           return;
1422         }
1423       }
1424     else //(mesh->GetType() == MED_STRUCTURED_MESH)
1425       {
1426       grid = vtkMedUnstructuredGrid::New();
1427       }
1428     grid->SetParentMesh(mesh);
1429     grid->SetComputeStep(cs);
1430     this->ReadGridInformation(grid);
1431     mesh->AddGridStep(grid);
1432     grid->Delete();
1433     grid->SetPreviousGrid(previousGrid);
1434     previousGrid = grid;
1435     }
1436 }
1437
1438 void vtkMedDriver30::ReadLocalizationInformation(vtkMedLocalization* loc)
1439 {
1440   FileOpen open(this);
1441
1442   med_int ngp;
1443   med_int spaceDimension;
1444   med_geometry_type type_geo;
1445   med_geometry_type sectiongeotype;
1446   med_int nsectionmeshcell;
1447
1448   char name[MED_NAME_SIZE+1] = "";
1449   char interpolationName[MED_NAME_SIZE+1] = "";
1450   char sectionName[MED_NAME_SIZE+1] = "";
1451
1452   if(MEDlocalizationInfo(
1453       this->FileId,
1454       loc->GetMedIterator(),
1455       name,
1456       &type_geo,
1457       &spaceDimension,
1458       &ngp,
1459       interpolationName,
1460       sectionName,
1461       &nsectionmeshcell,
1462       &sectiongeotype ) < 0)
1463     {
1464     vtkErrorMacro("Reading information on quadrature points definition : "
1465         << loc->GetMedIterator());
1466     }
1467
1468   loc->SetName(name);
1469   loc->SetInterpolationName(interpolationName);
1470   loc->SetSectionName(sectionName);
1471   loc->SetNumberOfQuadraturePoint(ngp);
1472   loc->SetGeometryType(type_geo);
1473   loc->SetSpaceDimension(spaceDimension);
1474   loc->SetNumberOfCellInSection(nsectionmeshcell);
1475   loc->SetSectionGeometryType(sectiongeotype);
1476
1477   med_float *localCoordinates = new med_float[loc->GetSizeOfPointLocalCoordinates()];
1478   med_float *pqLocalCoordinates = new med_float[loc->GetSizeOfQuadraturePointLocalCoordinates()];
1479   med_float *weights = new med_float[loc->GetSizeOfWeights()];
1480
1481   if(MEDlocalizationRd(FileId,
1482       loc->GetName(),
1483       MED_FULL_INTERLACE,
1484       localCoordinates,
1485       pqLocalCoordinates,
1486       weights) < 0)
1487     {
1488     vtkErrorMacro("MEDlocalizationRd : " << loc->GetName());
1489     }
1490
1491   vtkDoubleArray* lc = loc->GetPointLocalCoordinates();
1492   vtkDoubleArray *pqlc = loc->GetQuadraturePointLocalCoordinates();
1493   vtkDoubleArray *w = loc->GetWeights();
1494
1495   lc->SetNumberOfValues(loc->GetSizeOfPointLocalCoordinates());
1496   for(int i=0; i<loc->GetSizeOfPointLocalCoordinates(); i++)
1497     {
1498     lc->SetValue(i, localCoordinates[i]);
1499     }
1500
1501   pqlc->SetNumberOfValues(loc->GetSizeOfQuadraturePointLocalCoordinates());
1502   for(int i=0; i<loc->GetSizeOfQuadraturePointLocalCoordinates(); i++)
1503     {
1504     pqlc->SetValue(i, pqLocalCoordinates[i]);
1505     }
1506
1507   w->SetNumberOfValues(loc->GetSizeOfWeights());
1508   for(int i=0; i<loc->GetSizeOfWeights(); i++)
1509     {
1510     w->SetValue(i, weights[i]);
1511     }
1512 }
1513
1514 void vtkMedDriver30::ReadInterpolationInformation(vtkMedInterpolation* interp)
1515 {
1516
1517   med_geometry_type geotype;
1518   med_bool cellnode;
1519   med_int nbofbasisfunc;
1520   med_int nbofvariable;
1521   med_int maxdegree;
1522   med_int nmaxcoef;
1523
1524   char name[MED_NAME_SIZE+1] = "";
1525
1526   if(MEDinterpInfo (this->FileId,
1527                     interp->GetMedIterator(),
1528                     name,
1529                     &geotype, &cellnode, &nbofbasisfunc,
1530                     &nbofvariable, &maxdegree, &nmaxcoef) < 0)
1531     {
1532     vtkErrorMacro("MEDinterpInfo");
1533     return;
1534     }
1535
1536   interp->SetName(name);
1537   interp->SetGeometryType(geotype);
1538   interp->SetIsCellNode(cellnode);
1539   interp->SetNumberOfVariable(nbofvariable);
1540   interp->SetMaximumDegree(maxdegree);
1541   interp->SetMaximumNumberOfCoefficient(nmaxcoef);
1542   interp->AllocateNumberOfBasisFunction(nbofbasisfunc);
1543
1544   for(int basisid=0; basisid < nbofbasisfunc; basisid++)
1545     {
1546     vtkMedFraction* func = interp->GetBasisFunction(basisid);
1547     func->SetNumberOfVariable(nbofvariable);
1548
1549     med_int ncoef = MEDinterpBaseFunctionCoefSize (
1550         this->FileId,
1551         interp->GetName(),
1552         basisid+1);
1553     func->SetNumberOfCoefficients(ncoef);
1554
1555     if(ncoef <= 0 || nbofvariable <= 0)
1556       continue;
1557
1558     med_int *power = new med_int[nbofvariable * ncoef];
1559     med_float *coefficient = new med_float[ncoef];
1560
1561     if(MEDinterpBaseFunctionRd  (
1562         this->FileId,
1563         interp->GetName(),
1564         basisid+1,
1565         &ncoef,
1566         power,
1567         coefficient) < 0)
1568       {
1569       vtkErrorMacro("MEDinterpBaseFunctionRd");
1570       continue;
1571       }
1572     vtkDoubleArray* coeffs = func->GetCoefficients();
1573     for(int cid=0; cid < ncoef; cid++)
1574       {
1575       coeffs->SetValue(cid, coefficient[cid]);
1576       }
1577     vtkIntArray* powers = func->GetPowers();
1578     for(int pid=0; pid < ncoef*nbofvariable; pid++)
1579       {
1580       powers->SetValue(pid, power[pid]);
1581       }
1582
1583     delete[] power;
1584     delete[] coefficient;
1585     }
1586 }
1587
1588 void vtkMedDriver30::ReadSupportMeshInformation(
1589     vtkMedMesh* supportMesh)
1590 {
1591   FileOpen open(this);
1592
1593   char supportmeshname[MED_NAME_SIZE+1] = "";
1594   char description[MED_COMMENT_SIZE+1] = "";
1595   med_int spacedim;
1596   med_int meshdim;
1597   med_axis_type axistype;
1598   char axisname[3*MED_SNAME_SIZE+1] = "";
1599   char axisunit[3*MED_SNAME_SIZE+1] = "";
1600
1601   if(MEDsupportMeshInfo (this->FileId,
1602                          supportMesh->GetMedIterator(),
1603                          supportmeshname,
1604                          &spacedim,
1605                          &meshdim,
1606                          description,
1607                          &axistype,
1608                          axisname,
1609                          axisunit) < 0)
1610     {
1611     vtkErrorMacro("MEDsupportMeshInfo");
1612     }
1613
1614   supportMesh->SetName(supportmeshname);
1615   supportMesh->SetDescription(description);
1616   supportMesh->SetSpaceDimension(spacedim);
1617   supportMesh->SetMeshDimension(meshdim);
1618   supportMesh->SetAxisType(axistype);
1619   for(int dim = 0; dim < 3; dim++)
1620     {
1621     char axisname_dim[MED_SNAME_SIZE+1] = "";
1622     char axisunit_dim[MED_SNAME_SIZE+1] = "";
1623
1624     strncpy(axisname_dim, axisname+dim*MED_SNAME_SIZE, MED_SNAME_SIZE);
1625     strncpy(axisunit_dim, axisunit+dim*MED_SNAME_SIZE, MED_SNAME_SIZE);
1626
1627     supportMesh->GetAxisName()->SetValue(dim, axisname_dim);
1628     supportMesh->GetAxisUnit()->SetValue(dim, axisunit_dim);
1629     }
1630
1631   return;
1632 }
1633
1634 void vtkMedDriver30::LoadFamilyIds(vtkMedEntityArray* array)
1635 {
1636   if(array->IsFamilyIdsLoaded())
1637     return;
1638
1639   FileOpen open(this);
1640
1641   vtkMedGrid* grid = array->GetParentGrid();
1642
1643   vtkMedComputeStep cs = grid->GetComputeStep();
1644
1645   // first, find if the family ids are implicit or explicit
1646   med_bool changement, transformation;
1647
1648   med_int nfamid = MEDmeshnEntity(this->FileId,
1649                       grid->GetParentMesh()->GetName(),
1650                       cs.TimeIt,
1651                       cs.IterationIt,
1652                       array->GetEntity().EntityType,
1653                       array->GetEntity().GeometryType,
1654                       MED_FAMILY_NUMBER,
1655                       MED_NO_CMODE,
1656                       &changement,
1657                       &transformation);
1658
1659   if(nfamid == array->GetNumberOfEntity())
1660     {
1661
1662     vtkMedIntArray* famIds = vtkMedIntArray::New();
1663     array->SetFamilyIds(famIds);
1664     famIds->Delete();
1665
1666     famIds->SetNumberOfTuples(nfamid);
1667
1668     if ( MEDmeshEntityFamilyNumberRd(
1669             this->FileId,
1670             grid->GetParentMesh()->GetName(),
1671             cs.TimeIt,
1672             cs.IterationIt,
1673             array->GetEntity().EntityType,
1674             array->GetEntity().GeometryType,
1675             famIds->GetPointer(0) ) < 0)
1676       {
1677       vtkWarningMacro("Error loading the family ids of entity "
1678         << array->GetEntity().EntityType
1679         << " " << array->GetEntity().GeometryType
1680         << " on mesh " << grid->GetParentMesh()->GetName());
1681       array->SetFamilyIds(NULL);
1682       }
1683     }
1684   else
1685     {
1686     vtkDebugMacro("NumberOfEntity != Number of family ids");
1687     array->SetFamilyIds(NULL);
1688     }
1689
1690   array->ComputeFamilies();
1691 }
1692
1693 void vtkMedDriver30::LoadPointGlobalIds(vtkMedGrid* grid)
1694 {
1695   if(grid->IsPointGlobalIdsLoaded())
1696     return;
1697
1698   FileOpen open(this);
1699
1700   vtkMedIntArray* globalIds = vtkMedIntArray::New();
1701   grid->SetPointGlobalIds(globalIds);
1702   globalIds->Delete();
1703
1704   globalIds->SetNumberOfTuples(grid->GetNumberOfPoints());
1705
1706   if( MEDmeshEntityNumberRd (
1707         this->FileId,
1708         grid->GetParentMesh()->GetName(),
1709         grid->GetComputeStep().TimeIt,
1710         grid->GetComputeStep().IterationIt,
1711         MED_NODE,
1712         MED_NONE,
1713         globalIds->GetPointer(0) ) < 0)
1714     {
1715     grid->SetPointGlobalIds(NULL);
1716     }
1717 }
1718
1719 void vtkMedDriver30::LoadCoordinates(vtkMedGrid* grid)
1720 {
1721   if(grid->IsCoordinatesLoaded())
1722     return;
1723
1724   vtkMedRegularGrid* rgrid = vtkMedRegularGrid::SafeDownCast(grid);
1725   if(rgrid != NULL)
1726     {
1727     this->LoadRegularGridCoordinates(rgrid);
1728     return;
1729     }
1730
1731   vtkMedUnstructuredGrid* ugrid = vtkMedUnstructuredGrid::SafeDownCast(grid);
1732   vtkMedCurvilinearGrid* cgrid = vtkMedCurvilinearGrid::SafeDownCast(grid);
1733   if(ugrid == NULL && cgrid == NULL)
1734     {
1735     //TODO : deal with structured grids
1736     vtkWarningMacro("this kind of grid is not yet supported");
1737     return;
1738     }
1739
1740   if(grid->GetUsePreviousCoordinates())
1741     {
1742     vtkMedGrid* previousgrid = grid->GetPreviousGrid();
1743     if(previousgrid == NULL)
1744       {
1745       vtkErrorMacro("coordiantes have not changed, "
1746                     << "but there is no previous grid!");
1747       return;
1748       }
1749
1750     this->LoadCoordinates(previousgrid);
1751     if(ugrid != NULL)
1752       ugrid->SetCoordinates(vtkMedUnstructuredGrid::SafeDownCast(previousgrid)
1753                             ->GetCoordinates());
1754     if(cgrid != NULL)
1755       cgrid->SetCoordinates(vtkMedCurvilinearGrid::SafeDownCast(previousgrid)
1756                             ->GetCoordinates());
1757     }
1758   else
1759     {
1760
1761     FileOpen open(this);
1762
1763     vtkDataArray* coords = vtkMedUtilities::NewCoordArray();
1764     if(ugrid != NULL)
1765       ugrid->SetCoordinates(coords);
1766     if(cgrid != NULL)
1767       cgrid->SetCoordinates(coords);
1768     coords->Delete();
1769
1770     vtkMedComputeStep cs = grid->GetComputeStep();
1771
1772     coords->SetNumberOfComponents(grid->GetParentMesh()->GetSpaceDimension());
1773     coords->SetNumberOfTuples(grid->GetNumberOfPoints());
1774
1775     if ( MEDmeshNodeCoordinateRd( this->FileId,
1776                                   grid->GetParentMesh()->GetName(),
1777                                   cs.TimeIt,
1778                                   cs.IterationIt,
1779                                   MED_FULL_INTERLACE,
1780                                   (med_float*) coords->GetVoidPointer(0) ) < 0)
1781       {
1782       vtkErrorMacro("Load Coordinates for mesh "
1783           << grid->GetParentMesh()->GetName());
1784       }
1785     }
1786 }
1787
1788 void vtkMedDriver30::LoadProfile(vtkMedProfile* profile)
1789 {
1790   if(!profile || profile->IsLoaded())
1791     return;
1792
1793   FileOpen open(this);
1794
1795   vtkMedIntArray* indices = vtkMedIntArray::New();
1796   profile->SetIds(indices);
1797   indices->Delete();
1798
1799   indices->SetNumberOfTuples(profile->GetNumberOfElement());
1800
1801   char name[MED_NAME_SIZE+1] = "";
1802
1803   if( MEDprofileRd(this->FileId,
1804                    profile->GetName(),
1805                    indices->GetPointer(0) ) < 0)
1806     {
1807     vtkErrorMacro("Reading profile indices ");
1808     }
1809 }
1810
1811 void vtkMedDriver30::LoadConnectivity(vtkMedEntityArray* array)
1812 {
1813   if(array->IsConnectivityLoaded())
1814     return;
1815
1816   FileOpen open(this);
1817
1818   vtkMedGrid* grid = array->GetParentGrid();
1819
1820   grid = array->GetParentGrid();
1821
1822   const char* meshName = grid->GetParentMesh()->GetName();
1823
1824   vtkMedIntArray* conn = vtkMedIntArray::New();
1825   array->SetConnectivityArray(conn);
1826   conn->Delete();
1827
1828   vtkMedComputeStep cs = grid->GetComputeStep();
1829
1830   med_bool change;
1831   med_bool transformation;
1832
1833   if(array->GetEntity().GeometryType == MED_POLYGON)
1834     {
1835     // first check if we have points
1836     med_int connSize = MEDmeshnEntity(
1837                             this->FileId,
1838                             meshName,
1839                             cs.TimeIt,
1840                             cs.IterationIt,
1841                             array->GetEntity().EntityType,
1842                             MED_POLYGON,
1843                             MED_CONNECTIVITY,
1844                             array->GetConnectivity(),
1845                             &change,
1846                             &transformation);
1847
1848     if (connSize < 0)
1849       {
1850       vtkErrorMacro(<< "Error while reading polygons connectivity size."
1851                     << endl );
1852             return;
1853       }
1854
1855     conn->SetNumberOfTuples(connSize);
1856
1857     // How many polygon in the mesh in nodal connectivity mode
1858     // For the polygons, we get the size of array index
1859     med_int indexsize;
1860     if ((indexsize = MEDmeshnEntity(this->FileId,
1861                                     meshName,
1862                                     cs.TimeIt,
1863                                     cs.IterationIt,
1864                                     array->GetEntity().EntityType,
1865                                     MED_POLYGON,
1866                                     MED_INDEX_NODE,
1867                                     array->GetConnectivity(),
1868                                     &change,
1869                                     &transformation )) < 0)
1870       {
1871       vtkErrorMacro(<< "Error while reading polygons array index." << endl );
1872             return;
1873       }
1874
1875     vtkMedIntArray* index = vtkMedIntArray::New();
1876     array->SetFaceIndex(index);
1877     index->Delete();
1878
1879     index->SetNumberOfTuples( indexsize );
1880
1881     if ( MEDmeshPolygonRd(this->FileId,
1882                           meshName,
1883                           cs.TimeIt,
1884                           cs.IterationIt,
1885                           array->GetEntity().EntityType,
1886                           array->GetConnectivity(),
1887                           index->GetPointer(0),
1888                           conn->GetPointer(0) ) < 0)
1889       {
1890       vtkErrorMacro(<< "MEDmeshPolygonRd");
1891       return;
1892       }
1893     }
1894   else if(array->GetEntity().GeometryType == MED_POLYHEDRON)
1895     {
1896
1897     vtkIdType connSize = MEDmeshnEntity(this->FileId,
1898                                         meshName,
1899                                         grid->GetComputeStep().TimeIt,
1900                                         grid->GetComputeStep().IterationIt,
1901                                         array->GetEntity().EntityType,
1902                                         MED_POLYHEDRON,
1903                                         MED_CONNECTIVITY,
1904                                         array->GetConnectivity(),
1905                                         &change,
1906                                         &transformation);
1907     if (connSize < 0)
1908       {
1909       vtkErrorMacro(<< "Error while reading polyhedrons connectivity size."
1910                     << endl );
1911             return;
1912       }
1913
1914     conn->SetNumberOfTuples(connSize);
1915
1916     vtkMedIntArray* faceIndex = vtkMedIntArray::New();
1917     array->SetFaceIndex(faceIndex);
1918     faceIndex->Delete();
1919
1920     vtkMedIntArray* nodeIndex = vtkMedIntArray::New();
1921     array->SetNodeIndex(nodeIndex);
1922     nodeIndex->Delete();
1923
1924     vtkIdType np = array->GetNumberOfEntity() + 1;
1925     faceIndex->SetNumberOfTuples(np);
1926
1927     med_int nodeIndexSize;
1928
1929     if ((nodeIndexSize = MEDmeshnEntity(this->FileId,
1930                                         meshName,
1931                                         grid->GetComputeStep().TimeIt,
1932                                         grid->GetComputeStep().IterationIt,
1933                                         array->GetEntity().EntityType,
1934                                         MED_POLYHEDRON,
1935                                         MED_INDEX_NODE,
1936                                         array->GetConnectivity(),
1937                                         &change,
1938                                         &transformation )) < 0)
1939       {
1940       vtkErrorMacro(<< "Error while reading polygons array index." << endl );
1941             return;
1942       }
1943
1944     nodeIndex->SetNumberOfTuples(nodeIndexSize);
1945
1946     if (MEDmeshPolyhedronRd(this->FileId,
1947                             meshName,
1948                             cs.TimeIt,
1949                             cs.IterationIt,
1950                             array->GetEntity().EntityType,
1951                             array->GetConnectivity(),
1952                             faceIndex->GetPointer(0),
1953                             nodeIndex->GetPointer(0),
1954                             conn->GetPointer(0) ) < 0)
1955       {
1956       vtkErrorMacro("Error while reading connectivity of polyhedrons");
1957       return;
1958       }
1959
1960     }
1961   else
1962     {
1963     bool doReadConnectivity = true;
1964     if(array->GetConnectivity() == MED_NODAL)
1965       {
1966       if(array->GetEntity().EntityType == MED_STRUCT_ELEMENT)
1967         {
1968         std::cout << "  -- MED_STRUCT_ELEMENT --" << std::endl;
1969         if(array->GetStructElement() == NULL)
1970           {
1971           vtkErrorMacro("Entity type = MED_STRUCT_ELEMENT, but StructElement is not set!");
1972           return;
1973           }
1974         vtkIdType ntuple = array->GetNumberOfEntity()
1975                            * array->GetStructElement()->GetConnectivitySize();
1976
1977         conn->SetNumberOfTuples(ntuple);
1978         // particles are special : connectivity is not stored in the med file
1979         if(strcmp(array->GetStructElement()->GetName(), MED_PARTICLE_NAME) == 0 )
1980           {
1981           for(vtkIdType cellId = 0; cellId < ntuple; cellId++)
1982             {
1983             conn->SetValue(cellId, cellId+1);
1984             }
1985           doReadConnectivity = false;
1986           }
1987         }
1988       else
1989         {
1990         conn->SetNumberOfTuples(array->GetNumberOfEntity()
1991             * vtkMedUtilities::GetNumberOfPoint(
1992                 array->GetEntity().GeometryType));
1993         }
1994       }
1995     else
1996       {
1997       conn->SetNumberOfTuples(array->GetNumberOfEntity()
1998           * vtkMedUtilities::GetNumberOfSubEntity(
1999               array->GetEntity().GeometryType));
2000       }
2001
2002     if  (this->ParallelFileId == -1) // also (array->GetFilter() == NULL)
2003       {
2004       if ( (MEDmeshElementConnectivityRd(
2005             this->FileId,
2006             meshName,
2007             cs.TimeIt,
2008             cs.IterationIt,
2009             array->GetEntity().EntityType,
2010             array->GetEntity().GeometryType,
2011             array->GetConnectivity(),
2012             MED_FULL_INTERLACE,
2013             conn->GetPointer(0)) ) < 0)
2014         {
2015         vtkErrorMacro("Error while load connectivity of cells "
2016             << array->GetEntity().GeometryType);
2017         }
2018       }
2019     else
2020       {
2021       med_filter filter = MED_FILTER_INIT;
2022
2023       int    start;
2024       int    stride;
2025       int    count;
2026       int    blocksize;
2027       int    lastblocksize;
2028       array->GetFilter()->GetFilterSizes(start, stride, count,
2029                                    blocksize, lastblocksize );
2030
2031       med_int nbofconstituentpervalue = vtkMedUtilities::GetNumberOfNodes(
2032                                         array->GetEntity().GeometryType);
2033
2034       if ( MEDfilterBlockOfEntityCr( this->ParallelFileId,
2035               array->GetNumberOfEntity(),
2036             1, // one is for mesh elements, more than 1 is for fields
2037               nbofconstituentpervalue,
2038             MED_ALL_CONSTITUENT,
2039             MED_FULL_INTERLACE,
2040             MED_COMPACT_STMODE,
2041             MED_NO_PROFILE,
2042             (med_size)start,
2043             (med_size)stride,
2044             (med_size)count,
2045             (med_size)blocksize,
2046             (med_size)lastblocksize,
2047             &filter ) < 0 )
2048         {
2049         vtkErrorMacro("Filter creation ");
2050         }
2051
2052         if ( (MEDmeshElementConnectivityAdvancedRd(
2053               this->ParallelFileId,
2054               meshName,
2055               cs.TimeIt,
2056               cs.IterationIt,
2057               array->GetEntity().EntityType,
2058               array->GetEntity().GeometryType,
2059               array->GetConnectivity(),
2060               &filter,
2061               conn->GetPointer(0)) ) < 0)
2062           {
2063           vtkErrorMacro("Error while load connectivity of cells "
2064               << array->GetEntity().GeometryType);
2065           }
2066
2067       if ( MEDfilterClose( &filter ) < 0)
2068           {
2069         vtkErrorMacro("ERROR : filter closing ...");
2070           }
2071
2072       }
2073     }
2074 }
2075
2076 void vtkMedDriver30::LoadCellGlobalIds(vtkMedEntityArray* array)
2077 {
2078   if(array->IsGlobalIdsLoaded())
2079     return;
2080
2081   FileOpen open(this);
2082
2083   vtkMedIntArray* globalIds = vtkMedIntArray::New();
2084   array->SetGlobalIds(globalIds);
2085   globalIds->Delete();
2086
2087   globalIds->SetNumberOfTuples(array->GetNumberOfEntity());
2088
2089   vtkMedGrid* grid = array->GetParentGrid();
2090   vtkMedComputeStep cs = grid->GetComputeStep();
2091
2092   if( MEDmeshEntityNumberRd (
2093         this->FileId,
2094         grid->GetParentMesh()->GetName(),
2095         cs.TimeIt,
2096         cs.IterationIt,
2097         array->GetEntity().EntityType,
2098         array->GetEntity().GeometryType,
2099         globalIds->GetPointer(0) ) < 0)
2100     {
2101     array->SetGlobalIds(NULL);
2102     }
2103 }
2104
2105 void vtkMedDriver30::LoadField(vtkMedFieldOnProfile* fop, med_storage_mode mode)
2106 {
2107   FileOpen open(this);
2108
2109   vtkMedFieldOverEntity* fieldOverEntity = fop->GetParentFieldOverEntity();
2110   vtkMedFieldStep *step = fieldOverEntity->GetParentStep();
2111   vtkMedField* field = step->GetParentField();
2112   const vtkMedComputeStep& cs = step->GetComputeStep();
2113
2114   vtkDataArray* data = vtkMedUtilities::NewArray(field->GetDataType());
2115   fop->SetData(data);
2116   data->Delete();
2117
2118   med_int size;
2119   if(mode == MED_COMPACT_STMODE)
2120     {
2121     size = fop->GetNumberOfValues();
2122     }
2123   else
2124     {
2125     med_int profilesize;
2126     med_int nbofintegrationpoint;
2127     char profileName[MED_NAME_SIZE+1] = "";
2128     char localizationName[MED_NAME_SIZE+1] = "";
2129     size = MEDfieldnValueWithProfile(this->FileId,
2130                 field->GetName(),
2131                 cs.TimeIt,
2132                 cs.IterationIt,
2133                 fieldOverEntity->GetEntity().EntityType,
2134                 fieldOverEntity->GetEntity().GeometryType,
2135                 fop->GetMedIterator(),
2136                 MED_GLOBAL_STMODE,
2137                 profileName,
2138                 &profilesize,
2139                 localizationName,
2140                 &nbofintegrationpoint);
2141     }
2142
2143   if(fop->GetNumberOfIntegrationPoint() > 1)
2144     {
2145     size *= fop->GetNumberOfIntegrationPoint();
2146     }
2147
2148   data->SetNumberOfComponents(field->GetNumberOfComponent());
2149   data->SetNumberOfTuples(size);
2150   if  (this->ParallelFileId == -1)
2151     {
2152     if ( MEDfieldValueWithProfileRd(
2153           this->FileId,
2154           field->GetName(),
2155           cs.TimeIt,
2156           cs.IterationIt,
2157           fieldOverEntity->GetEntity().EntityType,
2158           fieldOverEntity->GetEntity().GeometryType,
2159           mode,
2160           fop->GetProfileName(),
2161           MED_FULL_INTERLACE,
2162           MED_ALL_CONSTITUENT,
2163           (unsigned char*) data->GetVoidPointer(0) ) < 0)
2164       {
2165       vtkErrorMacro("Error on MEDfieldValueWithProfileRd");
2166       }
2167     }
2168   else
2169     {
2170   if  (field->GetFieldType() == vtkMedField::CellField)
2171     {
2172     med_filter filter = MED_FILTER_INIT;
2173
2174     int    start;
2175     int    stride;
2176     int    count;
2177     int    blocksize;
2178     int    lastblocksize;
2179     fop->GetFilter()->GetFilterSizes(start, stride, count,
2180                                  blocksize, lastblocksize );
2181
2182     if ( MEDfilterBlockOfEntityCr( this->ParallelFileId,
2183         fop->GetNumberOfValues(),
2184           1, // one is for mesh elements, more than 1 is for fields
2185           field->GetNumberOfComponent(),
2186           MED_ALL_CONSTITUENT,
2187           MED_FULL_INTERLACE,
2188           MED_COMPACT_STMODE,
2189           MED_NO_PROFILE,
2190           (med_size)start,
2191           (med_size)stride,
2192           (med_size)count,
2193           (med_size)blocksize,
2194           (med_size)lastblocksize,
2195           &filter ) < 0 )
2196       {
2197       vtkErrorMacro("Filter creation ");
2198       }
2199
2200     if ( MEDfieldValueAdvancedRd(
2201             this->ParallelFileId,
2202             field->GetName(),
2203             cs.TimeIt,
2204             cs.IterationIt,
2205             fieldOverEntity->GetEntity().EntityType,
2206             fieldOverEntity->GetEntity().GeometryType,
2207             &filter,
2208             (unsigned char*) data->GetVoidPointer(0) ) < 0)
2209         {
2210         vtkErrorMacro("Error on MEDfieldValueAdvancedRd");
2211         }
2212
2213     if ( MEDfilterClose( &filter ) < 0)
2214         {
2215       vtkErrorMacro("ERROR : filter closing ...");
2216         }
2217       }
2218   else
2219     {//TODO : option utilisateur pour desactiver ou non les champs avec profile en //
2220     if ( MEDfieldValueWithProfileRd(
2221               this->FileId,
2222               field->GetName(),
2223               cs.TimeIt,
2224               cs.IterationIt,
2225               fieldOverEntity->GetEntity().EntityType,
2226               fieldOverEntity->GetEntity().GeometryType,
2227               mode,
2228               fop->GetProfileName(),
2229               MED_FULL_INTERLACE,
2230               MED_ALL_CONSTITUENT,
2231               (unsigned char*) data->GetVoidPointer(0) ) < 0)
2232           {
2233           vtkErrorMacro("Error on MEDfieldValueWithProfileRd");
2234           }
2235     }
2236     }
2237 }
2238
2239 void vtkMedDriver30::LoadVariableAttribute(vtkMedVariableAttribute* varatt,
2240                                            vtkMedEntityArray* array)
2241 {
2242   FileOpen open(this);
2243
2244   void  *value = NULL;
2245
2246   vtkAbstractArray* valuearray = array->GetVariableAttributeValue(varatt);
2247   // first test if this is already loaded
2248   if(valuearray != NULL && valuearray->GetNumberOfTuples() > 0)
2249     return;
2250
2251   if(valuearray == NULL)
2252     {
2253     valuearray = vtkMedUtilities::NewArray(varatt->GetAttributeType());
2254     array->SetVariableAttributeValues(varatt, valuearray);
2255     valuearray->Delete();
2256     }
2257
2258   valuearray->SetNumberOfComponents(varatt->GetNumberOfComponent());
2259   valuearray->SetNumberOfTuples(array->GetNumberOfEntity());
2260   valuearray->SetName(varatt->GetName());
2261
2262   vtkSmartPointer<vtkCharArray> chararray = vtkSmartPointer<vtkCharArray>::New();
2263
2264   if(varatt->GetAttributeType() != MED_ATT_NAME)
2265     {
2266     value = valuearray->GetVoidPointer(0);
2267     }
2268   else
2269     {
2270     chararray->SetNumberOfValues(varatt->GetNumberOfComponent() *
2271                                   array->GetNumberOfEntity() *
2272                                   MED_NAME_SIZE);
2273
2274     value = chararray->GetVoidPointer(0);
2275     }
2276
2277   vtkMedComputeStep cs = array->GetParentGrid()->GetComputeStep();
2278
2279   if(MEDmeshStructElementVarAttRd(
2280       this->FileId,
2281       array->GetParentGrid()->GetParentMesh()->GetName(),
2282       cs.TimeIt,
2283       cs.IterationIt,
2284       varatt->GetParentStructElement()->GetGeometryType(),
2285       varatt->GetName(),
2286       value) < 0)
2287     {
2288
2289     if(cs.IterationIt == MED_NO_IT && cs.TimeIt == MED_NO_DT && cs.TimeOrFrequency == MED_UNDEF_DT)
2290       {
2291       vtkErrorMacro("MEDmeshStructElementVarAttRd");
2292       return;
2293       }
2294     // try to see if I can reuse
2295     // the variable attributes of the NO_DT, NO_IT compute step
2296     vtkMedComputeStep nocs;
2297     nocs.IterationIt = MED_NO_IT;
2298     nocs.TimeIt = MED_NO_DT;
2299     nocs.TimeOrFrequency = MED_UNDEF_DT;
2300     vtkMedEntityArray* nocs_array =
2301         array->GetParentGrid()->GetParentMesh()->GetGridStep(nocs)->GetEntityArray(array->GetEntity());
2302     if(nocs_array == NULL)
2303       {
2304       nocs_array = array->GetParentGrid()->GetParentMesh()->GetGridStep(0)->GetEntityArray(array->GetEntity());
2305       }
2306
2307     if(nocs_array == NULL || nocs_array == array)
2308       {
2309       // try to force load the default compute step.
2310       if(MEDmeshStructElementVarAttRd(
2311           this->FileId,
2312           array->GetParentGrid()->GetParentMesh()->GetName(),
2313           nocs.TimeIt,
2314           nocs.IterationIt,
2315           varatt->GetParentStructElement()->GetGeometryType(),
2316           varatt->GetName(),
2317           value) < 0)
2318         {
2319         vtkErrorMacro("MEDmeshStructElementVarAttRd");
2320         return;
2321         }
2322       }
2323     else
2324       {
2325       this->LoadVariableAttribute(varatt, nocs_array);
2326       array->SetVariableAttributeValues(varatt, nocs_array->GetVariableAttributeValue(varatt));
2327       return;
2328       }
2329     }
2330
2331   // If I am here, it means that I read the values
2332   if(varatt->GetAttributeType() == MED_ATT_NAME)
2333     {
2334     char current_name[MED_NAME_SIZE+1] = "";
2335     vtkStringArray* sarray = vtkStringArray::SafeDownCast(valuearray);
2336     for(vtkIdType id = 0; id < varatt->GetNumberOfComponent() *
2337                        array->GetNumberOfEntity(); id++)
2338       {
2339       memset(current_name, '\0', MED_NAME_SIZE+1);
2340       strncpy(current_name, ((char*)value) + id*MED_NAME_SIZE, MED_NAME_SIZE);
2341       sarray->SetValue(id, current_name);
2342       }
2343     }
2344
2345   return;
2346 }
2347
2348 void vtkMedDriver30::PrintSelf(ostream& os, vtkIndent indent)
2349 {
2350   this->Superclass::PrintSelf(os, indent);
2351 }