src/MEDMEM/MEDMEM_VtkMedDriver.cxx

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