src/MEDMEM/MEDMEM_VtkMedDriver.cxx

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