MeshFormatReader.cxx

   1 // Copyright (C) 2021-2023  CEA, EDF
   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, or (at your option) any later version.
   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 "MeshFormatReader.hxx"
  21
  22 #include "MEDFileMesh.hxx"
  23 #include "MEDFileField.hxx"
  24 #include "MEDFileData.hxx"
  25 #include "MEDCouplingFieldDouble.hxx"
  26 #include "libmesh5.hxx"
  27 #include "MEDMESHConverterUtilities.hxx"
  28 #include <cstring>
  29 #include <fstream>
  30
  31 namespace MEDCoupling {
  32 MeshFormatReader::MeshFormatReader():_myMeshName("MESH")
  33 {}
  34 MeshFormatReader::MeshFormatReader(const std::string& meshFileName,
  35                                    const std::vector<std::string>& fieldFileName):_myFile(meshFileName),
  36                                    _myFieldFileNames(fieldFileName),
  37                                    _myCurrentFileId(0),
  38                                    _myMeshName("MESH")
  39 {}
  40
  41 MeshFormatReader::~MeshFormatReader()
  42 {}
  43
  44 MEDCoupling::MCAuto<MEDCoupling::MEDFileData> MeshFormatReader::loadInMedFileDS()
  45 {
  46     _myStatus = perform();
  47     if(_myStatus != MeshFormat::DRS_OK) return 0;
  48
  49     if ( !_uMesh->getName().c_str() || strlen( _uMesh->getName().c_str() ) == 0 )
  50         _uMesh->setName( _myMeshName );
  51     if ( _myFieldFileNames.size() ) performFields();
  52
  53
  54     MEDCoupling::MCAuto< MEDCoupling::MEDFileMeshes > meshes = MEDCoupling::MEDFileMeshes::New();
  55     _myMed = MEDCoupling::MEDFileData::New();
  56     meshes->pushMesh( _uMesh );
  57     _myMed->setMeshes( meshes );
  58
  59     if ( _fields ) _myMed->setFields( _fields );
  60
  61     return _myMed.retn();
  62 }
  63
  64 //================================================================================
  65 /*!
  66  * \brief Read a GMF file
  67  */
  68 //================================================================================
  69
  70 MeshFormat::Status MeshFormatReader::perform()
  71 {
  72     MeshFormat::Localizer loc;
  73
  74     MeshFormat::Status status = MeshFormat::DRS_OK;
  75
  76
  77     // open the file
  78      _reader = MeshFormat::MeshFormatParser();
  79     _myCurrentOpenFile = _myFile;
  80     _myCurrentFileId = _reader.GmfOpenMesh( _myFile.c_str(), GmfRead, &_version, &_dim );
  81     if ( !_myCurrentFileId )
  82     {
  83         if ( MeshFormat::isMeshExtensionCorrect( _myFile ))
  84         {
  85
  86             return addMessage( MeshFormat::Comment("Can't open for reading ") << _myFile,
  87                                /*fatal=*/true );
  88         }
  89         else
  90             return addMessage( MeshFormat::Comment("Not '.mesh' or '.meshb' extension of file ") << _myFile,
  91                                /*fatal=*/true );
  92     }
  93
  94
  95     // Read nodes
  96
  97
  98     MEDCoupling::DataArrayDouble* coordArray = MEDCoupling::DataArrayDouble::New();
  99     setNodes(coordArray);
 100
 101
 102     int nbEdges = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfEdges);
 103     int nbTria = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfTriangles);
 104     int nbQuad = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfQuadrilaterals);
 105     int nbTet = _reader.GmfStatKwd( _myCurrentFileId, MeshFormat::GmfTetrahedra );
 106     int nbPyr = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfPyramids);
 107     int nbHex = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfHexahedra);
 108     int nbPrism = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfPrisms);
 109
 110     _dim1NbEl = nbEdges;
 111     _dim2NbEl = nbTria + nbQuad;
 112     _dim3NbEl = nbTet + nbPyr + nbHex + nbPrism;
 113     bool okdim1 = (nbEdges > 0), okdim2= (_dim2NbEl > 0), okdim3= (_dim3NbEl > 0);
 114
 115     MEDCoupling::MEDCouplingUMesh* dimMesh1;
 116     MEDCoupling::MEDCouplingUMesh* dimMesh2;
 117     MEDCoupling::MEDCouplingUMesh* dimMesh3;
 118
 119     // dim 1
 120     if (okdim1 ) {
 121         dimMesh1 = MEDCoupling::MEDCouplingUMesh::New();
 122         dimMesh1->setCoords( coordArray );
 123         dimMesh1->allocateCells( nbEdges);
 124         dimMesh1->setMeshDimension( 1 );
 125     }
 126     // dim 2
 127     if (okdim2) {
 128         dimMesh2 = MEDCoupling::MEDCouplingUMesh::New();
 129         dimMesh2->setCoords( coordArray );
 130         dimMesh2->allocateCells(_dim2NbEl);
 131         dimMesh2->setMeshDimension( 2 );
 132     }
 133     // dim 3
 134     if (okdim3) {
 135         dimMesh3 = MEDCoupling::MEDCouplingUMesh::New();
 136         dimMesh3->setCoords( coordArray );
 137         dimMesh3->allocateCells(_dim3NbEl);
 138         dimMesh3->setMeshDimension( 3 );
 139     }
 140
 141     // Read elements
 142
 143     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 144
 145     /* Read edges */
 146
 147     if ( nbEdges)
 148     {
 149         setEdges(dimMesh1, nbEdges);
 150         dimMesh1->decrRef();
 151     }
 152
 153     /* Read triangles */
 154     if (nbTria)
 155     {
 156         setTriangles(dimMesh2, nbTria);
 157     }
 158
 159     /* Read quadrangles */
 160     if (nbQuad)
 161     {
 162         setQuadrangles( dimMesh2, nbQuad);
 163     }
 164
 165     if (okdim2 ) {
 166         dimMesh2->finishInsertingCells();
 167         _uMesh->setMeshAtLevel( 2 - _dim, dimMesh2 );
 168         dimMesh2->sortCellsInMEDFileFrmt();
 169         dimMesh2->decrRef();
 170     }
 171     /* Read terahedra */
 172     if (  nbTet )
 173     {
 174         setTetrahedras( dimMesh3, nbTet);
 175     }
 176
 177     /* Read pyramids */
 178     if ( nbPyr )
 179     {
 180         setPyramids( dimMesh3, nbPyr);
 181     }
 182
 183     /* Read hexahedra */
 184     if ( nbHex )
 185     {
 186         setHexahedras(dimMesh3, nbHex);
 187     }
 188
 189     /* Read prism */
 190     //~const int prismIDShift = myMesh->GetMeshInfo().NbElements();
 191     if ( nbPrism )
 192     {
 193         setPrisms(dimMesh3, nbPrism);
 194     }
 195
 196     if (okdim3 ) {
 197         dimMesh3->finishInsertingCells();
 198         _uMesh->setMeshAtLevel( 3 - _dim, dimMesh3 );
 199         dimMesh3->decrRef();
 200     }
 201
 202     buildFamilies();
 203     coordArray->decrRef();
 204     _reader.GmfCloseMesh(_myCurrentFileId);
 205     _myCurrentFileId = -1;
 206     _myCurrentOpenFile ="";
 207     return status;
 208 }
 209
 210 MeshFormat::Status MeshFormatReader::performFields()
 211 {
 212
 213     MeshFormat::Status status = MeshFormat::DRS_OK;
 214     _fields = MEDCoupling::MEDFileFields::New();
 215
 216
 217     const MeshFormat::GmfKwdCod meshFormatSol[1] = { MeshFormat::GmfSolAtVertices};  // ONLY VERTEX SOL FOR NOW
 218
 219     int dim, version;
 220     std::vector<std::string>::const_iterator fieldFileIt = _myFieldFileNames.begin();
 221
 222     for (; fieldFileIt !=_myFieldFileNames.end();  ++fieldFileIt)
 223     {
 224         _myCurrentOpenFile = *fieldFileIt;
 225         _myCurrentFileId = _reader.GmfOpenMesh( fieldFileIt->c_str(), GmfRead, &version, &dim );
 226         if ( !_myCurrentFileId )
 227         {
 228             if ( MeshFormat::isMeshExtensionCorrect( *fieldFileIt ))
 229             {
 230
 231                 return addMessage( MeshFormat::Comment("Can't open for reading ") << *fieldFileIt,
 232                                    /*fatal=*/true );
 233             }
 234             else
 235                 return addMessage( MeshFormat::Comment("Not '.sol' or '.solb' extension of file ") << _myFile,
 236                                    /*fatal=*/true );
 237         }
 238
 239         if (version != _version)
 240         {
 241
 242             addMessage( MeshFormat::Comment("Warning sol file version is different than mesh file version ") << *fieldFileIt,
 243                                /*fatal=*/false );
 244         }
 245         if (dim != _dim)
 246         {
 247
 248             return addMessage( MeshFormat::Comment("Error sol file must have same dimension as mesh file ") << *fieldFileIt,
 249                                /*fatal=*/true );
 250         }
 251
 252
 253         MeshFormat::GmfKwdCod kwd = meshFormatSol[0];
 254         int NmbSol, NmbTypes, NmbReals, TypesTab[ GmfMaxTyp ];
 255         NmbSol = _reader.GmfStatKwd( _myCurrentFileId, kwd, &NmbTypes, &NmbReals, TypesTab );
 256         if(NmbSol)
 257         {
 258
 259             for(int i=0; i<NmbTypes; i++)
 260             {
 261                 _reader.GmfGotoKwd(_myCurrentFileId, kwd);
 262                 switch(TypesTab[i])
 263                 {
 264                 case GmfSca:
 265                 {
 266
 267                     setFields(kwd, NmbSol, 1);
 268                     break;
 269                 }
 270                 case GmfVec:
 271                 {
 272
 273                     setFields(kwd, NmbSol, dim);
 274                     break;
 275                 }
 276                 case GmfSymMat :
 277                 {
 278
 279                     setFields(kwd, NmbSol, dim*(dim+1)/2 );
 280                     break;
 281                 }
 282                 case GmfMat:
 283                 {
 284
 285                     setFields(kwd, NmbSol, dim*dim);
 286                     break;
 287                 }
 288                 }
 289
 290             }
 291
 292
 293
 294         }
 295
 296
 297
 298         _reader.GmfCloseMesh(_myCurrentFileId);
 299         _myCurrentFileId = -1;
 300         _myCurrentOpenFile ="";
 301     }
 302     return status;
 303 }
 304
 305 void MeshFormatReader::setFields( MeshFormat::GmfKwdCod kwd, int nmbSol, int nbComp)
 306 {
 307
 308     MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> fieldValues =  MEDCoupling::DataArrayDouble::New();
 309     fieldValues->alloc(nmbSol, nbComp);
 310     double* values = fieldValues->getPointer();
 311
 312     int ref;
 313     double *val = new double[nbComp];
 314
 315     bool isOnAll = (_uMesh->getNumberOfNodes()== nmbSol);
 316
 317
 318     for(int i = 1; i<= nmbSol; i++)
 319     {
 320         callParserGetLin(kwd, val, nbComp, &ref);
 321
 322         std::copy(val, val+nbComp, values);
 323
 324         values+=nbComp;
 325     }
 326     values=fieldValues->getPointer(); // to delete
 327
 328     delete []val;
 329
 330     MEDCoupling::MCAuto<MEDCoupling::MEDCouplingFieldDouble> timeStamp;
 331     MEDCoupling::MCAuto<MEDCoupling::MEDFileFieldMultiTS> tsField =  MEDCoupling::MEDFileFieldMultiTS::New();
 332     MEDCoupling::MCAuto<MEDCoupling::MEDCouplingUMesh> dimMesh ;
 333     int  dimRel;
 334     MEDCoupling::TypeOfField typeOfField;
 335     setTypeOfFieldAndDimRel(kwd,  &typeOfField,  &dimRel );
 336
 337     timeStamp = MEDCoupling::MEDCouplingFieldDouble::New(typeOfField);
 338     dimMesh = _uMesh->getMeshAtLevel(dimRel);
 339
 340
 341     timeStamp->setMesh( dimMesh );
 342     std::string name = "Field_on_Vertex";
 343     timeStamp->setName(name);
 344     timeStamp->setArray(fieldValues);
 345     // set an order
 346     const int nbTS = tsField->getNumberOfTS();
 347     if ( nbTS > 0 )
 348         timeStamp->setOrder( nbTS );
 349
 350     // add the time-stamp
 351     timeStamp->checkConsistencyLight();
 352
 353     if (isOnAll)
 354     {
 355         tsField->appendFieldNoProfileSBT( timeStamp );
 356     }
 357
 358
 359     _fields->pushField( tsField);
 360
 361 }
 362
 363 void MeshFormatReader::setMeshName(const std::string& theMeshName)
 364 {
 365     _myMeshName = theMeshName;
 366 }
 367
 368 std::string MeshFormatReader::getMeshName() const
 369 {
 370     return _myMeshName;
 371 }
 372
 373
 374 void MeshFormatReader::setFile(const std::string& theFileName)
 375 {
 376     _myFile = theFileName;
 377 }
 378
 379 void MeshFormatReader::setFieldFileNames(const std::vector<std::string>& theFieldFileNames)
 380 {
 381     _myFieldFileNames = theFieldFileNames;
 382 }
 383 std::vector<std::string> MeshFormatReader::getFieldFileNames() const
 384 {
 385     return _myFieldFileNames;
 386 }
 387
 388 MeshFormat::Status MeshFormatReader::addMessage(const std::string& msg,
 389         const bool         isFatal/*=false*/)
 390 {
 391     if ( isFatal )
 392         _myErrorMessages.clear(); // warnings are useless if a fatal error encounters
 393
 394     _myErrorMessages.push_back( msg );
 395
 396     //~MESSAGE(msg);
 397 #ifdef _DEBUG_
 398     std::cout << msg << std::endl;
 399 #endif
 400     return ( _myStatus = isFatal ? MeshFormat::DRS_FAIL : MeshFormat::DRS_WARN_SKIP_ELEM );
 401 }
 402
 403 void MeshFormatReader::backward_shift(mcIdType* tab, int size)
 404 {
 405     for(int i = 0; i<size; i++) tab[i] = tab[i]-1;
 406 }
 407
 408
 409 MeshFormat::Status MeshFormatReader::setNodes( MEDCoupling::DataArrayDouble* coordArray)
 410 {
 411
 412     MeshFormat::Status status = MeshFormat::DRS_OK;
 413     int nbNodes = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfVertices);
 414     if ( nbNodes < 1 )
 415         return addMessage( "No nodes in the mesh", /*fatal=*/true );
 416
 417     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfVertices);
 418
 419     int ref;
 420     _uMesh = MEDCoupling::MEDFileUMesh::New();
 421     coordArray->alloc( nbNodes, _dim );
 422     double* coordPrt = coordArray->getPointer();
 423     std::vector<double> nCoords (_dim+2, 0.0) ;
 424
 425     double* coordPointer  = &nCoords[0];
 426
 427     if ( _version != GmfFloat )
 428     {
 429         double x, y, z;
 430
 431         for ( int i = 1; i <= nbNodes; ++i )
 432         {
 433             _dim == 2 ? _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfVertices, &x, &y, &ref) : \
 434             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfVertices, &x, &y, &z, &ref);
 435             nCoords[0] = x;
 436             nCoords[1] = y;
 437             nCoords[2] = z;
 438             std::copy(coordPointer, coordPointer+_dim, coordPrt);
 439             MeshFormatElement e(MeshFormat::GmfVertices, i-1);
 440             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 1);
 441             coordPrt += _dim;
 442
 443         }
 444
 445     }
 446     else
 447     {
 448         float x, y, z;
 449
 450         for ( int i = 1; i <= nbNodes; ++i )
 451         {
 452             _dim == 2 ? _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfVertices, &x, &y, &ref) :
 453             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfVertices, &x, &y, &z, &ref);
 454             nCoords[0] = x;
 455             nCoords[1] = y;
 456             nCoords[2] = z;
 457             std::copy(coordPointer, coordPointer+_dim, coordPrt);
 458             MeshFormatElement e(MeshFormat::GmfVertices, i-1);
 459             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 1);
 460             coordPrt += _dim;
 461         }
 462     }
 463     _uMesh->setCoords( coordArray );
 464
 465     return status;
 466
 467 }
 468
 469
 470 void MeshFormatReader::setEdges( MEDCoupling::MEDCouplingUMesh* dimMesh1, int nbEdges)
 471 {
 472
 473     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 474     int ref;
 475     // read extra vertices for quadratic edges
 476     std::vector<int> quadNodesAtEdges( nbEdges + 1, -1 );
 477     if ( int nbQuadEdges = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtEdges))
 478     {
 479         _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtEdges);
 480         for ( int i = 1; i <= nbQuadEdges; ++i )
 481         {
 482             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtEdges, &iN[0], &iN[1], &iN[2]);
 483             if ( iN[1] >= 1 )
 484                 quadNodesAtEdges[ iN[0] ] = iN[2];
 485         }
 486     }
 487     // create edges
 488     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfEdges);
 489     for ( int i = 1; i <= nbEdges; ++i )
 490     {
 491         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfEdges, &iN[0], &iN[1], &ref);
 492         const int midN = quadNodesAtEdges[ i ];
 493         if ( midN > 0 )
 494         {
 495
 496             mcIdType nodalConnPerCell[3] = {iN[0], iN[1], midN};
 497             backward_shift(nodalConnPerCell, 3);
 498             dimMesh1->insertNextCell(INTERP_KERNEL::NORM_SEG3, 3,nodalConnPerCell);
 499             MeshFormatElement e(MeshFormat::GmfEdges, i-1);
 500             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 1 -_dim);
 501         }
 502         else
 503         {
 504
 505             mcIdType nodalConnPerCell[2] = {iN[0], iN[1]};
 506             backward_shift(nodalConnPerCell, 2);
 507             dimMesh1->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2,nodalConnPerCell);
 508             MeshFormatElement e(MeshFormat::GmfEdges, i-1);
 509             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 1 -_dim);
 510         }
 511     }
 512
 513     dimMesh1->finishInsertingCells();
 514     dimMesh1->sortCellsInMEDFileFrmt();
 515     _uMesh->setMeshAtLevel( 1 - _dim, dimMesh1 );
 516 }
 517
 518 void MeshFormatReader::setTriangles(MEDCoupling::MEDCouplingUMesh* dimMesh2, int nbTria)
 519 {
 520     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 521     int ref;
 522     // read extra vertices for quadratic triangles
 523     std::vector< std::vector<int> > quadNodesAtTriangles( nbTria + 1 );
 524     if ( int nbQuadTria = _reader.GmfStatKwd(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtTriangles ))
 525     {
 526         _reader.GmfGotoKwd( _myCurrentFileId, MeshFormat::GmfExtraVerticesAtTriangles );
 527         for ( int i = 1; i <= nbQuadTria; ++i )
 528         {
 529             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtTriangles,
 530                              &iN[0], &iN[1], &iN[2], &iN[3], &iN[4],
 531                              &iN[5]); // iN[5] - preview TRIA7
 532             if ( iN[0] <= nbTria )
 533             {
 534                 std::vector<int>& nodes = quadNodesAtTriangles[ iN[0] ];
 535                 nodes.insert( nodes.end(), & iN[2], & iN[5+1] );
 536                 nodes.resize( iN[1] );
 537             }
 538         }
 539     }
 540
 541     // create triangles
 542
 543     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfTriangles);
 544
 545     for ( int i = 1; i <= nbTria; ++i )
 546     {
 547         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfTriangles, &iN[0], &iN[1], &iN[2], &ref);
 548         std::vector<int>& midN = quadNodesAtTriangles[ i ];
 549         if ( midN.size() >= 3 )
 550         {
 551
 552             mcIdType nodalConnPerCell[6] = {iN[0], iN[1], iN[2], midN[0], midN[1], midN[2]};
 553             backward_shift(nodalConnPerCell, 6);
 554             dimMesh2->insertNextCell(INTERP_KERNEL::NORM_TRI6, 6, nodalConnPerCell);
 555             MeshFormatElement e(MeshFormat::GmfTriangles, i-1);
 556             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 2 -_dim);
 557         }
 558         else
 559         {
 560             mcIdType nodalConnPerCell[3] = {iN[0], iN[1], iN[2]};
 561             backward_shift(nodalConnPerCell, 3);
 562             dimMesh2->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,nodalConnPerCell);
 563             MeshFormatElement e(MeshFormat::GmfTriangles, i-1);
 564             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 2 -_dim);
 565         }
 566         if ( !midN.empty() ) MeshFormat::FreeVector( midN );
 567     }
 568
 569
 570 }
 571
 572 void MeshFormatReader::setQuadrangles( MEDCoupling::MEDCouplingUMesh* dimMesh2, int nbQuad)
 573 {
 574     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 575     int ref;
 576     // read extra vertices for quadratic quadrangles
 577     std::vector< std::vector<int> > quadNodesAtQuadrilaterals( nbQuad + 1 );
 578     if ( int nbQuadQuad = _reader.GmfStatKwd( _myCurrentFileId, MeshFormat::GmfExtraVerticesAtQuadrilaterals ))
 579     {
 580         _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtQuadrilaterals);
 581         for ( int i = 1; i <= nbQuadQuad; ++i )
 582         {
 583             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtQuadrilaterals,
 584                              &iN[0], &iN[1], &iN[2], &iN[3], &iN[4], &iN[5], &iN[6]);
 585             if ( iN[0] <= nbQuad )
 586             {
 587                 std::vector<int>& nodes = quadNodesAtQuadrilaterals[ iN[0] ];
 588                 nodes.insert( nodes.end(), & iN[2], & iN[6+1] );
 589                 nodes.resize( iN[1] );
 590             }
 591         }
 592     }
 593     // create quadrangles
 594     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfQuadrilaterals);
 595     for ( int i = 1; i <= nbQuad; ++i )
 596     {
 597         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfQuadrilaterals, &iN[0], &iN[1], &iN[2], &iN[3], &ref);
 598         std::vector<int>& midN = quadNodesAtQuadrilaterals[ i ];
 599         if ( midN.size() == 8-4 ) // QUAD8
 600         {
 601             mcIdType nodalConnPerCell[8] = {iN[0], iN[1], iN[2], iN[3],
 602                                             midN[0], midN[1], midN[2], midN[3]
 603                                            };
 604             backward_shift(nodalConnPerCell, 8);
 605             dimMesh2->insertNextCell(INTERP_KERNEL::NORM_QUAD8,8, nodalConnPerCell);
 606             MeshFormatElement e(MeshFormat::GmfQuadrilaterals, i-1);
 607             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 2 -_dim);
 608         }
 609         else if ( midN.size() > 8-4 ) // QUAD9
 610         {
 611             mcIdType nodalConnPerCell[9] = {iN[0], iN[1], iN[2], iN[3],
 612                                             midN[0], midN[1], midN[2], midN[3], midN[4]
 613                                            };
 614             backward_shift(nodalConnPerCell, 9);
 615             dimMesh2->insertNextCell(INTERP_KERNEL::NORM_QUAD9,9, nodalConnPerCell);
 616             MeshFormatElement e(MeshFormat::GmfQuadrilaterals, i-1);
 617             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 2 -_dim);
 618         }
 619         else // QUAD4
 620         {
 621             mcIdType nodalConnPerCell[4] = {iN[0], iN[1], iN[2], iN[3]};
 622             backward_shift(nodalConnPerCell, 4);
 623             dimMesh2->insertNextCell(INTERP_KERNEL::NORM_QUAD4, 4, nodalConnPerCell);
 624             MeshFormatElement e(MeshFormat::GmfQuadrilaterals, i-1);
 625             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 2 -_dim);
 626         }
 627         if ( !midN.empty() ) MeshFormat::FreeVector( midN );
 628     }
 629
 630 }
 631
 632 void MeshFormatReader::setTetrahedras( MEDCoupling::MEDCouplingUMesh* dimMesh3, int nbTet)
 633 {
 634     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 635     int ref;
 636     // read extra vertices for quadratic tetrahedra
 637     std::vector< std::vector<int> > quadNodesAtTetrahedra( nbTet + 1 );
 638     if ( int nbQuadTetra = _reader.GmfStatKwd( _myCurrentFileId, MeshFormat::GmfExtraVerticesAtTetrahedra ))
 639     {
 640         _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtTetrahedra);
 641         for ( int i = 1; i <= nbQuadTetra; ++i )
 642         {
 643             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtTetrahedra,
 644                              &iN[0], &iN[1], &iN[2], &iN[3], &iN[4], &iN[5], &iN[6], &iN[7]);
 645             if ( iN[0] <= nbTet )
 646             {
 647                 std::vector<int>& nodes = quadNodesAtTetrahedra[ iN[0] ];
 648                 nodes.insert( nodes.end(), & iN[2], & iN[7+1] );
 649                 nodes.resize( iN[1] );
 650             }
 651         }
 652     }
 653     // create tetrahedra
 654     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfTetrahedra);
 655     for ( int i = 1; i <= nbTet; ++i )
 656     {
 657         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfTetrahedra, &iN[0], &iN[1], &iN[2], &iN[3], &ref);
 658         std::vector<int>& midN = quadNodesAtTetrahedra[ i ];
 659         if ( midN.size() >= 10-4 ) // TETRA10
 660         {
 661             mcIdType nodalConnPerCell[10] = {iN[0], iN[2], iN[1], iN[3],
 662                                              midN[2], midN[1], midN[0], midN[3], midN[5], midN[4]
 663                                             };
 664             backward_shift(nodalConnPerCell, 10);
 665             dimMesh3->insertNextCell(INTERP_KERNEL::NORM_TETRA10, 10,nodalConnPerCell);
 666             MeshFormatElement e(MeshFormat::GmfTetrahedra, i-1);
 667             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 668         }
 669         else // TETRA4
 670         {
 671             mcIdType nodalConnPerCell[4] = {iN[0], iN[2], iN[1], iN[3]};
 672             backward_shift(nodalConnPerCell, 4);
 673             dimMesh3->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4, nodalConnPerCell);
 674             MeshFormatElement e(MeshFormat::GmfTetrahedra, i-1);
 675             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 676         }
 677         if ( !midN.empty() ) MeshFormat::FreeVector( midN );
 678     }
 679
 680 }
 681
 682
 683 void MeshFormatReader::setPyramids( MEDCoupling::MEDCouplingUMesh* dimMesh3, int nbPyr)
 684 {
 685     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 686     int ref;
 687     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfPyramids);
 688     for ( int i = 1; i <= nbPyr; ++i )
 689     {
 690         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfPyramids, &iN[0], &iN[1], &iN[2], &iN[3], &iN[4], &ref);
 691         mcIdType nodalConnPerCell[5] = {iN[3], iN[2], iN[1], iN[0], iN[4]};
 692         backward_shift(nodalConnPerCell, 5);
 693         dimMesh3->insertNextCell(INTERP_KERNEL::NORM_PYRA5, 5,nodalConnPerCell);
 694         MeshFormatElement e(MeshFormat::GmfPyramids, i-1);
 695         _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 696     }
 697
 698 }
 699
 700
 701 void MeshFormatReader::setHexahedras( MEDCoupling::MEDCouplingUMesh* dimMesh3, int nbHex)
 702 {
 703     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 704     int ref;
 705     // read extra vertices for quadratic hexahedra
 706     std::vector< std::vector<int> > quadNodesAtHexahedra( nbHex + 1 );
 707     if ( int nbQuadHexa = _reader.GmfStatKwd( _myCurrentFileId, MeshFormat::GmfExtraVerticesAtHexahedra ))
 708     {
 709         _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtHexahedra);
 710         for ( int i = 1; i <= nbQuadHexa; ++i )
 711         {
 712             _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfExtraVerticesAtHexahedra, &iN[0], &iN[1], // Hexa Id, Nb extra vertices
 713                              &iN[2], &iN[3], &iN[4], &iN[5],
 714                              &iN[6], &iN[7], &iN[8], &iN[9],
 715                              &iN[10], &iN[11], &iN[12], &iN[13], // HEXA20
 716                              &iN[14],
 717                              &iN[15], &iN[16], &iN[17], &iN[18],
 718                              &iN[19],
 719                              &iN[20]);                          // HEXA27
 720             if ( iN[0] <= nbHex )
 721             {
 722                 std::vector<int>& nodes = quadNodesAtHexahedra[ iN[0] ];
 723                 nodes.insert( nodes.end(), & iN[2], & iN[20+1] );
 724                 nodes.resize( iN[1] );
 725             }
 726         }
 727     }
 728     // create hexhedra
 729
 730
 731     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfHexahedra);
 732     for ( int i = 1; i <= nbHex; ++i )
 733     {
 734         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfHexahedra, &iN[0], &iN[1], &iN[2], &iN[3],
 735                          &iN[4], &iN[5], &iN[6], &iN[7], &ref);
 736         std::vector<int>& midN = quadNodesAtHexahedra[ i ];
 737         if ( midN.size() == 20-8 ) // HEXA20
 738         {
 739             mcIdType nodalConnPerCell[20] = {iN[0], iN[3], iN[2], iN[1],
 740                                              iN[4], iN[7], iN[6], iN[5],
 741                                              midN[3], midN[2], midN[1], midN[0],
 742                                              midN[7], midN[6], midN[5], midN[4],
 743                                              midN[8], midN[11], midN[10], midN[9]
 744                                             };
 745             backward_shift(nodalConnPerCell, 20);
 746             dimMesh3->insertNextCell(INTERP_KERNEL::NORM_HEXA20,20, nodalConnPerCell);
 747             MeshFormatElement e(MeshFormat::GmfHexahedra, i-1);
 748             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 749
 750
 751         }
 752         else if ( midN.size() >= 27-8 ) // HEXA27
 753         {
 754             mcIdType nodalConnPerCell[27] = {iN[0], iN[3], iN[2], iN[1],
 755                                              iN[4], iN[7], iN[6], iN[5],
 756                                              midN[3], midN[2], midN[1], midN[0],
 757                                              midN[7], midN[6], midN[5], midN[4],
 758                                              midN[8], midN[11], midN[10], midN[9],
 759                                              midN[12],
 760                                              midN[16], midN[15], midN[14], midN[13],
 761                                              midN[17],
 762                                              midN[18]
 763                                             };
 764             backward_shift(nodalConnPerCell, 27);
 765             dimMesh3->insertNextCell(INTERP_KERNEL::NORM_HEXA27,27, nodalConnPerCell);
 766             MeshFormatElement e(MeshFormat::GmfHexahedra, i-1);
 767             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 768
 769         }
 770         else // HEXA8
 771         {
 772             mcIdType nodalConnPerCell[8] = {iN[0], iN[3], iN[2], iN[1],
 773                                             iN[4], iN[7], iN[6], iN[5]
 774                                            };
 775             backward_shift(nodalConnPerCell, 8);
 776             dimMesh3->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8, nodalConnPerCell);
 777             MeshFormatElement e(MeshFormat::GmfHexahedra, i-1);
 778             _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 779
 780         }
 781         if ( !midN.empty() ) MeshFormat::FreeVector( midN );
 782     }
 783
 784
 785 }
 786
 787
 788
 789 void MeshFormatReader::setPrisms(MEDCoupling::MEDCouplingUMesh* dimMesh3, int nbPrism)
 790 {
 791     int iN[28]; // 28 - nb nodes in HEX27 (+ 1 for safety :)
 792     int ref;
 793     _reader.GmfGotoKwd(_myCurrentFileId, MeshFormat::GmfPrisms);
 794     for ( int i = 1; i <= nbPrism; ++i )
 795     {
 796         _reader.GmfGetLin(_myCurrentFileId, MeshFormat::GmfPrisms, &iN[0], &iN[1], &iN[2], &iN[3], &iN[4], &iN[5], &ref);
 797         mcIdType nodalConnPerCell[8] = {iN[0], iN[2], iN[1], iN[3], iN[5], iN[4]};
 798         backward_shift(nodalConnPerCell, 8);
 799         dimMesh3->insertNextCell(INTERP_KERNEL::NORM_PENTA6, 6,nodalConnPerCell);
 800         MeshFormatElement e(MeshFormat::GmfPrisms, i-1);
 801         _fams.insert(std::pair <int, MeshFormatElement> (ref, e), 3 -_dim);
 802
 803     }
 804
 805
 806 }
 807
 808 void MeshFormatReader::callParserGetLin( MeshFormat::GmfKwdCod kwd,  double* val, int valSize, int* ref)
 809 {
 810     switch(valSize)
 811     {
 812     case 1:
 813     {
 814         _reader.GmfGetLin(_myCurrentFileId, kwd, &val[0], ref);
 815         break;
 816     }
 817     case 2:
 818     {
 819         _reader.GmfGetLin(_myCurrentFileId, kwd, &val[0], &val[1], ref);
 820         break;
 821     }
 822     case 3:
 823     {
 824         _reader.GmfGetLin(_myCurrentFileId, kwd, &val[0], &val[1], &val[2], ref);
 825         break;
 826     }
 827     case 4:
 828     {
 829         _reader.GmfGetLin(_myCurrentFileId, kwd, &val[0], &val[1], &val[2], &val[3], ref);
 830         break;
 831     }
 832     case 6:
 833     {
 834         _reader.GmfGetLin(_myCurrentFileId, kwd, &val[0], &val[1], &val[2], &val[3], &val[4], &val[5], ref);
 835         break;
 836     }
 837     case 9:
 838     {
 839         _reader.GmfGetLin(_myCurrentFileId, kwd, &val[0], &val[1], &val[2], &val[3], &val[4], &val[5], &val[6], &val[7], &val[8], ref);
 840         break;
 841     }
 842     }
 843 }
 844
 845 void MeshFormatReader::setTypeOfFieldAndDimRel(MeshFormat::GmfKwdCod kwd, MEDCoupling::TypeOfField* typeOfField, int* dimRel )
 846 {
 847     switch (kwd)
 848     {
 849     case MeshFormat::GmfSolAtVertices :
 850     {
 851         *typeOfField = MEDCoupling::ON_NODES;
 852         *dimRel = 1 ;
 853         break;
 854     }
 855     case MeshFormat::GmfSolAtEdges :
 856     {
 857         *typeOfField = MEDCoupling::ON_CELLS;
 858         *dimRel = 1 - _dim ;
 859         break;
 860     }
 861     case MeshFormat::GmfSolAtTriangles :
 862     case MeshFormat::GmfSolAtQuadrilaterals :
 863     {
 864         *typeOfField = MEDCoupling::ON_CELLS;
 865         *dimRel = 2 - _dim;
 866         break;
 867     }
 868     case MeshFormat::GmfSolAtTetrahedra :
 869     case MeshFormat::GmfSolAtPrisms :
 870     case MeshFormat::GmfSolAtHexahedra :
 871     {
 872         *typeOfField = MEDCoupling::ON_CELLS;
 873         *dimRel = 3 - _dim;
 874         break;
 875     }
 876     }
 877 }
 878
 879
 880 INTERP_KERNEL::NormalizedCellType MeshFormatReader::toMedType(MeshFormat::GmfKwdCod kwd)
 881 {
 882     INTERP_KERNEL::NormalizedCellType type;
 883     //~switch (kwd)
 884     //~{
 885     //~case MeshFormat::GmfEdges :
 886     //~{
 887
 888     //~type = INTERP_KERNEL::NORM_SEG2;
 889     //~break;
 890     //~}
 891     //~case MeshFormat::GmfTriangles :
 892     //~{
 893     //~type = INTERP_KERNEL::NORM_TRI3;
 894     //~break;
 895     //~}
 896     //~case MeshFormat::GmfQuadrilaterals :
 897     //~{
 898     //~type = INTERP_KERNEL::NORM_QUAD;
 899     //~break;
 900     //~}
 901     //~case MeshFormat::GmfSolAtTetrahedra :
 902     //~case MeshFormat::GmfSolAtPrisms :
 903     //~case MeshFormat::GmfSolAtHexahedra :
 904     //~{
 905     //~*typeOfField = MEDCoupling::ON_CELLS;
 906     //~*dimRel = 3 - _dim;
 907     //~break;
 908     //~}
 909     //~}
 910     return type;
 911 }
 912
 913
 914 void MeshFormatReader::buildFamilies()
 915 {
 916     buildNodesFamilies();
 917     buildCellsFamilies();
 918 }
 919
 920 void MeshFormatReader::buildCellsFamilies()
 921 {
 922     std::vector<int> levs =  _uMesh->getNonEmptyLevels();
 923     for (size_t iDim = 0; iDim<levs.size(); iDim++ )
 924     {
 925         int dimRelMax = levs[iDim];
 926         std::map <int, std::vector<MeshFormatElement>* > famDim = _fams.getMapAtLevel(dimRelMax);
 927         std::map <int, std::vector<MeshFormatElement>* >::const_iterator _meshFormatFamsIt = famDim.begin();
 928         std::vector< const MEDCoupling::DataArrayIdType* > fams;
 929         MEDCoupling::DataArrayIdType* cellIds = MEDCoupling::DataArrayIdType::New();
 930         cellIds->alloc(_uMesh->getSizeAtLevel(dimRelMax), 1);
 931         cellIds->fillWithZero();
 932
 933         for(; _meshFormatFamsIt!= famDim.end(); ++_meshFormatFamsIt)
 934         {
 935             const int famId = _meshFormatFamsIt->first;
 936             std::string famName ="FromMeshGemsFormatAttributFamily_"+std::to_string(famId);
 937             std::vector <MeshFormatElement>* cellsInFam = _meshFormatFamsIt->second;
 938             if (!famId) continue;
 939             std::vector <MeshFormatElement>::iterator cellsInFamIt = cellsInFam->begin();
 940
 941             _uMesh->addFamily(famName, famId);
 942             for ( ; cellsInFamIt !=cellsInFam->end(); ++cellsInFamIt)
 943             {
 944                 cellIds->setIJ(cellsInFamIt->_id, 0, famId);
 945             }
 946
 947         }
 948         _uMesh->setFamilyFieldArr(dimRelMax, cellIds->deepCopy());
 949         cellIds->decrRef();
 950
 951
 952     }
 953 }
 954
 955 void MeshFormatReader::buildNodesFamilies()
 956 {
 957   std::vector<int> levs =  _uMesh->getNonEmptyLevels();
 958      int dimRelMax = 1;
 959   std::map <int, std::vector<MeshFormatElement>* > famDim = _fams.getMapAtLevel(dimRelMax);
 960   std::map <int, std::vector<MeshFormatElement>* >::const_iterator _meshFormatFamsIt = famDim.begin();
 961   std::vector< const MEDCoupling::DataArrayIdType* > fams;
 962   MEDCoupling::DataArrayIdType* cellIds = MEDCoupling::DataArrayIdType::New();
 963   cellIds->alloc(_uMesh->getSizeAtLevel(dimRelMax), 1);
 964   cellIds->fillWithZero();
 965
 966   for(; _meshFormatFamsIt!= famDim.end(); ++_meshFormatFamsIt)
 967   {
 968     const int famId = _meshFormatFamsIt->first;
 969     if (!famId) continue;
 970     bool thisIsACellFamily = false;
 971
 972     for (size_t iDim = 0; iDim<levs.size(); iDim++ )
 973     {
 974       int dimMesh = levs[iDim];
 975       std::map <int, std::vector<MeshFormatElement>* > famDimAtLevel = _fams.getMapAtLevel(dimMesh);
 976       std::map <int, std::vector<MeshFormatElement>* >::iterator famDimAtLevelId = famDimAtLevel.find(famId);
 977       if (famDimAtLevelId != famDimAtLevel.end())
 978       {
 979         thisIsACellFamily = true;
 980         break;
 981       }
 982
 983     }
 984
 985     if (thisIsACellFamily) continue;
 986     std::string famName ="FromMeshGemsFormatAttributFamily_"+std::to_string(famId);
 987     std::vector <MeshFormatElement>* cellsInFam = _meshFormatFamsIt->second;
 988     std::vector <MeshFormatElement>::iterator cellsInFamIt = cellsInFam->begin();
 989
 990     _uMesh->addFamily(famName, famId);
 991     for ( ; cellsInFamIt !=cellsInFam->end(); ++cellsInFamIt)
 992     {
 993       cellIds->setIJ(cellsInFamIt->_id, 0, famId);
 994     }
 995
 996   }
 997   _uMesh->setFamilyFieldArr(dimRelMax, cellIds->deepCopy());
 998   cellIds->decrRef();
 999 }
1000 }