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