src/MEDMEM/MEDMEM_VtkMedDriver.cxx

   1 #include "MEDMEM_VtkMedDriver.hxx"
   2
   3 #include <sstream>
   4
   5 #include "MEDMEM_define.hxx"
   6 #include "MEDMEM_Med.hxx"
   7 #include "MEDMEM_Field.hxx"
   8 #include "MEDMEM_Support.hxx"
   9 #include "MEDMEM_Mesh.hxx"
  10 #include "MEDMEM_CellModel.hxx"
  11
  12 using namespace std;
  13 using namespace MEDMEM;
  14 using namespace MED_EN;
  15
  16 VTK_MED_DRIVER::VTK_MED_DRIVER(): GENDRIVER(),
  17                                   _ptrMed((MED * const)MED_NULL)
  18 {
  19   _vtkFile = new ofstream();
  20   // What about _id in Gendriver ?
  21   // _driverType ???
  22 }
  23
  24
  25 VTK_MED_DRIVER::VTK_MED_DRIVER(const string & fileName,  MED * const ptrMed):
  26   GENDRIVER(fileName,MED_EN::MED_RDWR), _ptrMed(ptrMed)
  27 {
  28   _ptrMed->addDriver(*this); // OU RECUPERER L'ID.
  29   _vtkFile = new ofstream();
  30   // What about _id in Gendriver ?
  31   // _driverType ???
  32 }
  33
  34 VTK_MED_DRIVER::VTK_MED_DRIVER(const VTK_MED_DRIVER & driver):
  35   GENDRIVER(driver),
  36   _ptrMed(driver._ptrMed)
  37 {
  38   _ptrMed->addDriver(*this); // OU RECUPERER L'ID.
  39   _vtkFile = new ofstream();
  40   // What about _id in Gendriver ?
  41   // _driverType ???
  42 }
  43
  44 VTK_MED_DRIVER::~VTK_MED_DRIVER()
  45 {
  46   const char * LOC ="VTK_MED_DRIVER::~VTK_MED_DRIVER()";
  47   BEGIN_OF(LOC);
  48
  49   close();
  50   delete _vtkFile ;
  51
  52   END_OF(LOC);
  53 }
  54
  55 GENDRIVER * VTK_MED_DRIVER::copy() const
  56 {
  57   return new VTK_MED_DRIVER(*this) ;
  58 }
  59
  60 //REM :  As t'on besoin du champ _status :  _vtkFile <-> _status  ?  Oui
  61
  62 void VTK_MED_DRIVER::openConst() const {
  63
  64   const char * LOC ="VTK_MED_DRIVER::open() : ";
  65   BEGIN_OF(LOC);
  66
  67   if ( _fileName == "" )
  68     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
  69                                      << "_fileName is |\"\"|, please set a correct fileName before calling open()"
  70                                      )
  71                           );
  72
  73   if (!(*_vtkFile).is_open())
  74     (*_vtkFile).open(_fileName.c_str()) ;
  75 //    if (*_vtkFile)
  76 //      _status = MED_OPENED ;
  77 //    else
  78   if (!(*_vtkFile))
  79     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Could not open file "
  80                                      << _fileName)
  81                           );
  82   END_OF(LOC);
  83 }
  84
  85 void VTK_MED_DRIVER::open() {
  86   openConst() ;
  87 }
  88
  89 void VTK_MED_DRIVER::closeConst() const {
  90
  91   const char * LOC = "VTK_MED_DRIVER::close() : ";
  92   BEGIN_OF(LOC);
  93
  94   (*_vtkFile).close();
  95
  96 //    if (*_vtkFile)
  97 //      _status = MED_CLOSED ;
  98 //    else
  99   if (!(*_vtkFile))
 100     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Could not close file "
 101                                      << _fileName)
 102                           );
 103   END_OF(LOC);
 104 }
 105
 106 void VTK_MED_DRIVER::close() {
 107   closeConst() ;
 108 }
 109
 110
 111 void VTK_MED_DRIVER::write() const {
 112
 113   const char * LOC = "VTK_MED_DRIVER::write() : ";
 114   BEGIN_OF(LOC);
 115
 116   // Well we must open vtk file first, because there are
 117   // no other driver than MED for VTK that do it !
 118   openConst() ;
 119
 120   // could we put more than one Mesh ?????
 121   (*_vtkFile) << "# vtk DataFile Version 2.0" << endl
 122            << "maillage from MedMemory"  << endl ;
 123   // only ASCII for the moment (binary came later :-)
 124   (*_vtkFile) << "ASCII" << endl ;
 125
 126   int NumberOfMeshes = _ptrMed->getNumberOfMeshes() ;
 127   deque<string> MeshNames = _ptrMed->getMeshNames() ;
 128   //deque<string>::const_iterator  currentMesh ; !! UNUSED VARIABLE !!
 129   // In fact, we must take care of all supports
 130   // We restrict Field on all nodes or cells
 131
 132   int NumberOfFields = _ptrMed->getNumberOfFields() ;
 133   deque<string> FieldNames = _ptrMed->getFieldNames() ;
 134   //deque<string>::const_iterator  currentField ; !! UNUSED VARIABLE !!
 135
 136   //  for ( currentMesh=MeshName.begin();currentMesh != MeshName.end(); currentMesh++) {
 137   for (int i=0; i<NumberOfMeshes; i++) {
 138     MESH * myMesh = _ptrMed->getMesh(MeshNames[i]) ;
 139     writeMesh(myMesh) ;
 140     // get all field which values are on this mesh => revoir api de Med !!!
 141     // first : field on node
 142     // fields is on all node !
 143     (*_vtkFile) << "POINT_DATA " << myMesh->getNumberOfNodes() << endl ;
 144     for (int j=0; j<NumberOfFields; j++) {
 145       deque<DT_IT_> timeStep = _ptrMed->getFieldIteration(FieldNames[j]) ;
 146       deque<DT_IT_>::const_iterator currentTimeStep ;
 147       for ( currentTimeStep=timeStep.begin(); currentTimeStep!=timeStep.end(); currentTimeStep++) {
 148         int dt = (*currentTimeStep).dt ;
 149         int it = (*currentTimeStep).it ;
 150         FIELD_ * myField = _ptrMed->getField(FieldNames[j],dt,it) ;
 151         if( MeshNames[i] == myField->getSupport()->getMesh()->getName() ) {
 152           // rigth in all case : better compare pointeur ?
 153           if (MED_NODE == myField->getSupport()->getEntity())
 154             if (myField->getSupport()->isOnAllElements()) {
 155               ostringstream name ;
 156               name << myField->getName() << "_" << dt << "_" << it ;
 157               writeField(myField,name.str()) ;
 158             } else
 159               MESSAGE(LOC << "Could not write field "<<myField->getName()<<" which is not on all nodes !");
 160         }
 161       }
 162     }
 163
 164     (*_vtkFile) << "CELL_DATA " << myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS) << endl ;
 165     // second : field on cell
 166     for (int j=0; j<NumberOfFields; j++) {
 167       deque<DT_IT_> timeStep = _ptrMed->getFieldIteration(FieldNames[j]) ;
 168       deque<DT_IT_>::const_iterator currentTimeStep ;
 169       for ( currentTimeStep=timeStep.begin(); currentTimeStep!=timeStep.end(); currentTimeStep++) {
 170         int dt ;
 171         int it ;
 172         FIELD_ * myField = _ptrMed->getField(FieldNames[j],dt,it) ;
 173         if( MeshNames[i] == myField->getSupport()->getMesh()->getName() ) {
 174           // rigth in all case : better compare pointeur ?
 175           if (MED_CELL == myField->getSupport()->getEntity())
 176             if (myField->getSupport()->isOnAllElements()) {
 177               ostringstream name ;
 178               name << myField->getName() << "_" << dt << "_" << it ;
 179               writeField(myField,name.str()) ;
 180             } else
 181               MESSAGE(LOC << "Could not write field "<<myField->getName()<<" which is not on all cells !");
 182         }
 183       }
 184     }
 185
 186   }
 187
 188   // Well we must close vtk file first, because there are
 189   // no other driver than MED for VTK that do it !
 190   //  closeConst() ;
 191
 192   END_OF(LOC);
 193 }
 194
 195 void VTK_MED_DRIVER::writeMesh(MESH * myMesh) const {
 196
 197   const char * LOC = "VTK_MED_DRIVER::writeMesh() : ";
 198   BEGIN_OF(LOC);
 199
 200   (*_vtkFile) << "DATASET UNSTRUCTURED_GRID" << endl ;
 201   // put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
 202   int SpaceDimension = myMesh->getSpaceDimension() ;
 203   int NumberOfNodes = myMesh->getNumberOfNodes() ;
 204   (*_vtkFile) << "POINTS " << NumberOfNodes << " float" << endl ;
 205   const double *coordinate = myMesh->getCoordinates(MED_FULL_INTERLACE) ;
 206   for (int i=0;i<NumberOfNodes;i++) {
 207     for (int j=0;j<SpaceDimension;j++)
 208       (*_vtkFile) << coordinate[i*SpaceDimension+j] << " " ;
 209     if (SpaceDimension==1)
 210       (*_vtkFile) << "0 0" ;
 211     if (SpaceDimension==2)
 212       (*_vtkFile) << "0" ;
 213     (*_vtkFile) << endl ;
 214   }
 215
 216   // we put connectivity
 217   // how many cells and how many value in connectivity :
 218   int cells_types_count = myMesh->getNumberOfTypes(MED_CELL) ;
 219   //  int * cells_count = myMesh->get_cells_count() ;
 220   //  int cells_sum = cells_count[cells_types_count] ;
 221   int cells_sum = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS) ;
 222   const CELLMODEL * cells_type = myMesh->getCellsTypes(MED_CELL) ;
 223   //  int connectivity_sum = 0 ;
 224
 225   //const int * connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS) ; !! UNUSED VARIABLE !!
 226   const int * connectivityIndex = myMesh->getConnectivityIndex(MED_NODAL,MED_CELL) ;
 227
 228   int connectivity_sum =  connectivityIndex[cells_sum]-1 ;
 229
 230   (*_vtkFile) << "CELLS " << cells_sum << " " << connectivity_sum+cells_sum << endl ;
 231   // we put connectivity
 232   for (int i=0;i<cells_types_count;i++) {
 233     int *filter = (int*) NULL ; // index in vtk connectivity
 234     switch (cells_type[i].getType())
 235       {
 236       case MED_POINT1  : {
 237         filter = new int[1] ;
 238         filter[0] = 0 ;
 239         break ;
 240       }
 241       case MED_SEG2    : {
 242         filter = new int[2] ;
 243         filter[0] = 0 ;
 244         filter[1] = 1 ;
 245         break ;
 246       }
 247       case MED_SEG3    : {
 248         break ;
 249       }
 250       case MED_TRIA3   : {
 251         filter = new int[3] ;
 252         filter[0] = 0 ;
 253         filter[1] = 1 ;
 254         filter[2] = 2 ;
 255         break ;
 256       }
 257       case MED_QUAD4   : {
 258         filter = new int[4] ;
 259         filter[0] = 0 ;
 260         filter[1] = 1 ;
 261         filter[2] = 2 ;
 262         filter[3] = 3 ;
 263         break ;
 264       }
 265       case MED_TRIA6   : {
 266         break ;
 267       }
 268       case MED_QUAD8   : {
 269         break ;
 270       }
 271       case MED_TETRA4  : {
 272         filter = new int[4] ;
 273         filter[0] = 0 ;
 274         filter[1] = 1 ;
 275         filter[2] = 3 ;  // 3td element in med are 4th in vtk (array begin at 0 !)
 276         filter[3] = 2 ;  // 4th element in med are 3rd in vtk (array begin at 0 !)
 277         break ;
 278       }
 279       case MED_PYRA5   : {
 280         filter = new int[5] ;
 281         filter[0] = 0 ;
 282         filter[1] = 3 ;  // 2nd element in med are 4th in vtk (array begin at 0 !)
 283         filter[2] = 2 ;
 284         filter[3] = 1 ;  // 4th element in med are 2nd in vtk (array begin at 0 !)
 285         filter[4] = 4 ;
 286         break ;
 287       }
 288       case MED_PENTA6  : {
 289         filter = new int[6] ;
 290         filter[0] = 0 ;
 291         filter[1] = 1 ;
 292         filter[2] = 2 ;
 293         filter[3] = 3 ;
 294         filter[4] = 4 ;
 295         filter[5] = 5 ;
 296         break ;
 297       }
 298       case MED_HEXA8   : {
 299         filter = new int[8] ;
 300         filter[0] = 0 ;
 301         filter[1] = 3 ;
 302         filter[2] = 2 ;
 303         filter[3] = 1 ;
 304         filter[4] = 4 ;
 305         filter[5] = 7 ;
 306         filter[6] = 6 ;
 307         filter[7] = 5 ;
 308         break ;
 309       }
 310       case MED_TETRA10 : {
 311         break ;
 312       }
 313       case MED_PYRA13  : {
 314         break ;
 315       }
 316       case MED_PENTA15 : {
 317         break ;
 318       }
 319       case MED_HEXA20  : {
 320         break ;
 321       }
 322       default : {
 323         break ;
 324       }
 325       }
 326     if (filter==NULL)
 327       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<": MED element type not supported yet : " << cells_type[i].getName() ) ) ;
 328     int nodes_cell = cells_type[i].getNumberOfNodes();
 329     int numberOfCell = myMesh->getNumberOfElements(MED_CELL,cells_type[i].getType()) ;
 330     const int * connectivityArray = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,cells_type[i].getType());
 331     for (int j=0;j<numberOfCell;j++) {
 332       (*_vtkFile) << nodes_cell << " " ;
 333       for (int k=0;k<nodes_cell;k++)
 334         (*_vtkFile) << connectivityArray[j*nodes_cell+filter[k]] - 1 << " " ;
 335       (*_vtkFile) << endl ;
 336     }
 337     if (filter != NULL)
 338       delete[] filter ;
 339   }
 340   (*_vtkFile) << endl ;
 341   // we put cells type
 342   (*_vtkFile) << "CELL_TYPES " << cells_sum << endl ;
 343   for (int i=0;i<cells_types_count;i++) {
 344     int vtkType = 0 ;
 345     switch (cells_type[i].getType())
 346       {
 347       case MED_POINT1  : {
 348         vtkType = 1 ;
 349         break ;
 350       }
 351       case MED_SEG2    : {
 352         vtkType = 3 ;
 353         break ;
 354       }
 355       case MED_SEG3    : {
 356         vtkType = 0 ;
 357         break ;
 358       }
 359       case MED_TRIA3   : {
 360         vtkType = 5 ;
 361         break ;
 362       }
 363       case MED_QUAD4   : {
 364         vtkType = 9 ;
 365         break ;
 366       }
 367       case MED_TRIA6   : {
 368         vtkType = 0 ;
 369         break ;
 370       }
 371       case MED_QUAD8   : {
 372         vtkType = 0 ;
 373         break ;
 374       }
 375       case MED_TETRA4  : {
 376         vtkType = 10 ;
 377         break ;
 378       }
 379       case MED_PYRA5   : {
 380         vtkType = 14 ;
 381         break ;
 382       }
 383       case MED_PENTA6  : {
 384         vtkType = 13 ;
 385         break ;
 386       }
 387       case MED_HEXA8   : {
 388         vtkType = 12 ;
 389         break ;
 390       }
 391       case MED_TETRA10 : {
 392         vtkType = 0 ;
 393         break ;
 394       }
 395       case MED_PYRA13  : {
 396         vtkType = 0 ;
 397         break ;
 398       }
 399       case MED_PENTA15 : {
 400         vtkType = 0 ;
 401         break ;
 402       }
 403       case MED_HEXA20  : {
 404         vtkType = 0 ;
 405         break ;
 406       }
 407       default : {
 408         vtkType = 0 ;
 409         break ;
 410       }
 411       }
 412     if (vtkType == 0)
 413       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<": MED element type not supported yet : " << cells_type[i].getType() ) ) ;
 414     int numberOfCell = myMesh->getNumberOfElements(MED_CELL,cells_type[i].getType()) ;
 415     for (int j=0;j<numberOfCell;j++)
 416       (*_vtkFile) << vtkType << endl ;
 417   }
 418
 419   // add a constant field on all node to test !
 420   //    _vtkFile << "POINT_DATA " << NumberOfNodes << endl ;
 421   //    _vtkFile << "SCALARS example_scalaire float 1" << endl ;
 422   //    _vtkFile << "LOOKUP_TABLE default" << endl ;
 423   //    for (int i=0;i<NumberOfNodes;i++)
 424   //      _vtkFile << i << endl ;
 425
 426   return ;
 427
 428
 429   END_OF(LOC);
 430 }
 431
 432 void VTK_MED_DRIVER::writeField(FIELD_ * myField,string name) const {
 433
 434   const char * LOC = "VTK_MED_DRIVER::writeField() : ";
 435   BEGIN_OF(LOC);
 436
 437   int NomberOfValue = myField->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS) ;
 438   int NomberOfComponents =  myField->getNumberOfComponents() ;
 439
 440   med_type_champ type = myField->getValueType() ;
 441   SCRUTE(name);
 442   SCRUTE(type);
 443   switch (type)
 444     {
 445     case MED_INT32 : {
 446       MESSAGE("MED_INT32");
 447       if (NomberOfComponents==3) {
 448         (*_vtkFile) << "VECTORS " << name << " int" << endl ;
 449       } else if (NomberOfComponents<=4) {
 450         (*_vtkFile) << "SCALARS " << name << " int " << NomberOfComponents << endl ;
 451         (*_vtkFile) << "LOOKUP_TABLE default" << endl ;
 452       } else {
 453         MESSAGE(LOC << "Could not write field "<<myField->getName()<<" there are more than 4 components !");
 454         return ;
 455       }
 456
 457       //const int * value = ((FIELD<int>*)myField)->getValue(MED_NO_INTERLACE) ;
 458       const int * value = ((FIELD<int>*)myField)->getValue(MED_NO_INTERLACE) ;
 459       for (int i=0; i<NomberOfValue; i++) {
 460         for(int j=0; j<NomberOfComponents; j++)
 461           (*_vtkFile) << value[j*NomberOfValue+i] << " " ;
 462         (*_vtkFile) << endl ;
 463       }
 464       break ;
 465     }
 466     case MED_REEL64 : {
 467       MESSAGE("MED_REEL64");
 468       if (NomberOfComponents==3) {
 469         (*_vtkFile) << "VECTORS " << name << " float" << endl ;
 470       } else if (NomberOfComponents<=4) {
 471         (*_vtkFile) << "SCALARS " << name << " float " << NomberOfComponents << endl ;
 472         (*_vtkFile) << "LOOKUP_TABLE default" << endl ;
 473       } else {
 474         MESSAGE(LOC << "Could not write field "<<myField->getName()<<" there are more than 4 components !");
 475         return ;
 476       }
 477       const double * value = ((FIELD<double>*)myField)->getValue(MED_NO_INTERLACE) ;
 478       for (int i=0; i<NomberOfValue; i++) {
 479         for(int j=0; j<NomberOfComponents; j++)
 480           (*_vtkFile) << value[j*NomberOfValue+i] << " " ;
 481         (*_vtkFile) << endl ;
 482       }
 483       break ;
 484     }
 485     default : {
 486       MESSAGE(LOC << "Could not write field "<<name<<" the type is not int or double !");
 487     }
 488     }
 489
 490   END_OF(LOC);
 491 }
 492
 493 void VTK_MED_DRIVER::writeSupport(SUPPORT * mySupport) const {
 494   const char * LOC = "VTK_MED_DRIVER::writeSupport(SUPPORT *) : " ;
 495   BEGIN_OF(LOC) ;
 496   MESSAGE(LOC << "Not yet implemented, acting on the object " << *mySupport);
 497   END_OF(LOC) ;
 498 }