src/DriverMED/DriverMED_W_Field.cxx

   1 // Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 // File      : DriverMED_W_Field.cxx
  24 // Created   : Thu Feb 27 17:45:00 2014
  25 // Author    : eap
  26
  27 #include "DriverMED_W_Field.h"
  28
  29 #include "DriverMED.hxx"
  30 #include "DriverMED_W_SMESHDS_Mesh.h"
  31 #include "MED_Factory.hxx"
  32 #include "MED_TFile.hxx"
  33 #include "MED_Utilities.hxx"
  34 #include "MED_Wrapper.hxx"
  35 #include "SMDS_IteratorOnIterators.hxx"
  36 #include "SMDS_MeshElement.hxx"
  37 #include "SMDS_SetIterator.hxx"
  38 #include "SMESHDS_Mesh.hxx"
  39 #include "SMESH_TypeDefs.hxx"
  40
  41 //================================================================================
  42 /*!
  43  * \brief Constructor
  44  */
  45 //================================================================================
  46
  47 DriverMED_W_Field::DriverMED_W_Field():
  48   //_medFileID( -1 ),
  49   _elemType( SMDSAbs_All ),
  50   _dt( -1 ),
  51   _it( -1 )
  52 {
  53 }
  54
  55 //================================================================================
  56 /*!
  57  * \brief Sets basic data
  58  *  \param [in] mesh - supporting mesh
  59  *  \param [in] fieldName - name of a field
  60  *  \param [in] type - type of supporting elements
  61  *  \param [in] nbComps - number of components
  62  *  \param [in] isIntData - type of data: double or integer
  63  */
  64 //================================================================================
  65
  66 bool DriverMED_W_Field::Set(SMESHDS_Mesh *      mesh,
  67                             const std::string & fieldName,
  68                             SMDSAbs_ElementType type,
  69                             const int           nbComps,
  70                             const bool          isIntData)
  71 {
  72   _fieldName = fieldName;
  73   _compNames.resize( nbComps, "" );
  74
  75   if ( type == SMDSAbs_All )
  76   {
  77     if ( mesh->NbVolumes() > 0 )
  78       type = SMDSAbs_Volume;
  79     else if ( mesh->NbFaces() > 0 )
  80       type = SMDSAbs_Face;
  81     else if ( mesh->NbEdges() > 0 )
  82       type = SMDSAbs_Edge;
  83     else
  84       type = SMDSAbs_Node;
  85   }
  86   if ( myMesh != mesh )
  87   {
  88     _nbElemsByGeom.clear();
  89     for ( int iG = 0; iG < SMDSEntity_Last; ++iG )
  90       _elemsByGeom[iG].clear();
  91     SetMesh( mesh );
  92   }
  93
  94   // find out "MED order" of elements - sort elements by geom type
  95   smIdType nbElems;
  96   if ( _nbElemsByGeom.empty() || _elemType != type )
  97   {
  98     _elemType = type;
  99     _nbElemsByGeom.resize( 1, std::make_pair( SMDSEntity_Last, 0 ));
 100
 101     // count nb of elems of each geometry
 102     SMDSAbs_EntityType geom = SMDSEntity_0D;
 103     for ( ; geom < SMDSEntity_Last; SMESHUtils::Increment( geom ))
 104     {
 105       SMDSAbs_ElementType t = SMDS_MeshCell::ElemType( geom );
 106       if ( t != _elemType ) continue;
 107
 108       nbElems = mesh->GetMeshInfo().NbElements( geom );
 109       if ( nbElems < 1 ) continue;
 110
 111       _nbElemsByGeom.push_back( std::make_pair( geom, nbElems + _nbElemsByGeom.back().second ));
 112     }
 113     // add nodes of missing 0D elements on VERTEXes
 114     if ( _addODOnVertices && _elemType == SMDSAbs_0DElement )
 115     {
 116       std::vector< const SMDS_MeshElement* >& nodes = _elemsByGeom[SMDSEntity_Node];
 117       if ( nodes.empty() )
 118         DriverMED_W_SMESHDS_Mesh::getNodesOfMissing0DOnVert( myMesh, nodes );
 119       if ( !nodes.empty() )
 120       {
 121         if ( _nbElemsByGeom.size() == 1 )
 122           _nbElemsByGeom.push_back( std::make_pair( SMDSEntity_0D, 0));
 123         _nbElemsByGeom.push_back( std::make_pair( SMDSEntity_Node,
 124                                                   nodes.size() + _nbElemsByGeom.back().second ));
 125       }
 126     }
 127
 128     // sort elements by their geometry
 129     int iGeoType, nbGeomTypes = _nbElemsByGeom.size() - 1;
 130     if ( nbGeomTypes > 1 )
 131     {
 132       for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 133       {
 134         iGeoType = _nbElemsByGeom[iG].first;
 135         nbElems  = _nbElemsByGeom[iG].second - _nbElemsByGeom[iG-1].second;
 136         _elemsByGeom[ iGeoType ].reserve( nbElems );
 137       }
 138       iGeoType = _nbElemsByGeom[1].first; // for missing 0D
 139       if ( _elemsByGeom[ iGeoType ].empty() )
 140       {
 141         nbElems = mesh->GetMeshInfo().NbElements( _elemType );
 142         SMDS_ElemIteratorPtr eIt = mesh->elementsIterator( _elemType );
 143         for ( int iE = 0; iE < nbElems && eIt->more(); ++iE )
 144         {
 145           const SMDS_MeshElement* e = eIt->next();
 146           _elemsByGeom[ e->GetEntityType() ].push_back( e );
 147         }
 148       }
 149     }
 150   }
 151   _intValues.clear();
 152   _dblValues.clear();
 153
 154   // allocate memory for values
 155   nbElems = _nbElemsByGeom.empty() ? 0 : _nbElemsByGeom.back().second;
 156   if ( isIntData )
 157     _intValues.reserve( nbElems * nbComps );
 158   else
 159     _dblValues.reserve( nbElems * nbComps );
 160
 161   return nbElems && nbComps;
 162 }
 163
 164 //================================================================================
 165 /*!
 166  * \brief Set a name of a component countered from zero
 167  */
 168 //================================================================================
 169
 170 void DriverMED_W_Field::SetCompName(const int iComp, const char* name)
 171 {
 172   if ( (int)_compNames.size() <= iComp )
 173     _compNames.resize( iComp + 1 );
 174   _compNames[ iComp ] = name;
 175 }
 176
 177 //================================================================================
 178 /*!
 179  * \brief Sets numdt and numit field features. Call this fun before AddValue()!
 180  */
 181 //================================================================================
 182
 183 void DriverMED_W_Field::SetDtIt(const int dt, const int it)
 184 {
 185   _dt = dt;
 186   _it = it;
 187   _intValues.clear();
 188   _dblValues.clear();
 189 }
 190
 191 //================================================================================
 192 /*!
 193  * \brief Adds a float field value
 194  */
 195 //================================================================================
 196
 197 void DriverMED_W_Field::AddValue( double val )
 198 {
 199   _dblValues.push_back( val );
 200 }
 201
 202 //================================================================================
 203 /*!
 204  * \brief Adds an integer field value
 205  */
 206 //================================================================================
 207
 208 void DriverMED_W_Field::AddValue( int    val )
 209 {
 210   _intValues.push_back( val );
 211 }
 212
 213 //================================================================================
 214 /*!
 215  * Returns elements in the order they are written in MED file
 216  */
 217 //================================================================================
 218
 219 SMDS_ElemIteratorPtr DriverMED_W_Field::GetOrderedElems()
 220 {
 221   if ( _nbElemsByGeom.size() < 2 )
 222     return SMDS_ElemIteratorPtr();
 223
 224   if ( _nbElemsByGeom.size() == 2 )
 225     // sole geom type of elements
 226     return myMesh->elementsIterator( _elemType );
 227
 228   std::vector< SMDS_ElemIteratorPtr > iterVec( _nbElemsByGeom.size()-1 );
 229   for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 230   {
 231     int iGeoType = _nbElemsByGeom[ iG ].first;
 232     iterVec[ iG-1 ] = SMDS_ElemIteratorPtr
 233       ( new SMDS_ElementVectorIterator( _elemsByGeom[ iGeoType ].begin(),
 234                                         _elemsByGeom[ iGeoType ].end() ));
 235   }
 236   typedef SMDS_IteratorOnIterators
 237     < const SMDS_MeshElement *, std::vector< SMDS_ElemIteratorPtr > > TItIterator;
 238   return SMDS_ElemIteratorPtr( new TItIterator( iterVec ));
 239 }
 240
 241 Driver_Mesh::Status DriverMED_W_Field::PerformInternal(MED::PWrapper& medFile)
 242 {
 243   if ( myMeshId < 0 && myMeshName.empty() )
 244     return addMessage("Mesh in file not specified", /*isFatal=*/true );
 245   if ( _nbElemsByGeom.size() < 2 )
 246     return addMessage("No values to write", /*isFatal=*/false );
 247   if ( !myMesh )
 248     return addMessage("Supporting mesh not set", /*isFatal=*/true );
 249
 250   MED::PMeshInfo meshInfo;
 251   if ( myMeshId > 0 )
 252   {
 253     meshInfo = medFile->GetPMeshInfo( myMeshId );
 254   }
 255   else
 256   {
 257     // look for a mesh by name
 258     int aNbMeshes = medFile->GetNbMeshes();
 259     for ( int iMesh = aNbMeshes; iMesh > 0 && myMeshId < 1; --iMesh )
 260     {
 261       meshInfo = medFile->GetPMeshInfo( iMesh );
 262       if ( !meshInfo || meshInfo->GetName() != myMeshName )
 263         meshInfo.reset();
 264       else
 265         myMeshId = iMesh;
 266     }
 267   }
 268   if (( !meshInfo ) ||
 269       ( !myMeshName.empty() && meshInfo->GetName() != myMeshName ))
 270   {
 271     myMeshId = -1;
 272     return addMessage("DriverMED_W_Field: Specified mesh not found in the file", /*isFatal=*/true );
 273   }
 274
 275   // create a field
 276   MED::ETypeChamp  dataType = _dblValues.empty() ? MED::eINT : MED::eFLOAT64;
 277   MED::PFieldInfo fieldInfo = medFile->CrFieldInfo( meshInfo,
 278                                                     _compNames.size(),
 279                                                     dataType );
 280   fieldInfo->SetName( _fieldName );
 281   for ( size_t iC = 0; iC < _compNames.size(); ++iC )
 282   {
 283     fieldInfo->SetCompName( iC, _compNames[ iC ]);
 284     fieldInfo->SetUnitName( iC, "");
 285   }
 286   if ( _compNames.size() > 1 )
 287   {
 288     for ( size_t i = 0; i < fieldInfo->myCompNames.size()-1; ++i )
 289       if ( !fieldInfo->myCompNames[i] )
 290         fieldInfo->myCompNames[i] = ' ';
 291   }
 292   medFile->SetFieldInfo( fieldInfo );
 293
 294   // specific treatment of added 0D elements
 295   if ( _nbElemsByGeom.size()   == 3 &&
 296        _nbElemsByGeom[1].first == SMDSEntity_0D )
 297   {
 298     _nbElemsByGeom[1].second += _nbElemsByGeom[2].second;
 299     _nbElemsByGeom.resize( 2 );
 300   }
 301
 302   // create a time stamp
 303   MED::TGeom2Size type2nb;
 304   for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 305   {
 306     SMDSAbs_EntityType    smdsType = _nbElemsByGeom[iG].first;
 307     MED::EGeometrieElement medType = (MED::EGeometrieElement) DriverMED::GetMedGeoType( smdsType );
 308     int                    nbElems = _nbElemsByGeom[iG].second - _nbElemsByGeom[iG-1].second;
 309     type2nb.insert( std::make_pair( medType, nbElems ));
 310   }
 311
 312   MED::EEntiteMaillage       entity = ( _elemType == SMDSAbs_Node ? MED::eNOEUD : MED::eMAILLE );
 313   MED::PTimeStampInfo timeStampInfo = medFile->CrTimeStampInfo( fieldInfo, entity, type2nb );
 314   timeStampInfo->myNumDt  = _dt;
 315   timeStampInfo->myNumOrd = _it;
 316
 317   MED::PTimeStampValueBase  timeStampVal = medFile->CrTimeStampValue( timeStampInfo, dataType );
 318   MED::PFloatTimeStampValue timeStampFltVal = timeStampVal;
 319   MED::PIntTimeStampValue   timeStampIntVal = timeStampVal;
 320
 321   // set values
 322   int iVal = 0;
 323   MED::TFloat* ptrDbl = 0;
 324   MED::TInt*   ptrInt = 0;
 325   for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 326   {
 327     SMDSAbs_EntityType    smdsType = _nbElemsByGeom[iG].first;
 328     MED::EGeometrieElement medType = (MED::EGeometrieElement) DriverMED::GetMedGeoType( smdsType );
 329     int nbElems = ( _nbElemsByGeom[iG].second - _nbElemsByGeom[iG-1].second ) * _compNames.size();
 330     if ( dataType == MED::eFLOAT64 )
 331     {
 332       ptrDbl = timeStampFltVal->GetMeshValue( medType ).GetPointer();
 333       for ( int i = 0; i < nbElems; ++i, ++iVal )
 334         ptrDbl[ i ] = _dblValues[ iVal ];
 335     }
 336     else
 337     {
 338       ptrInt = timeStampIntVal->GetMeshValue( medType ).GetPointer();
 339       for ( int i = 0; i < nbElems; ++i, ++iVal )
 340         ptrInt[ i ] = _intValues[ iVal ];
 341     }
 342   }
 343
 344   // write
 345   medFile->SetTimeStampValue( timeStampVal );
 346
 347   _dblValues.clear();
 348   _intValues.clear();
 349
 350   return DRS_OK;
 351 }
 352
 353 /*!
 354  * Writes a field to the file
 355  */
 356 Driver_Mesh::Status DriverMED_W_Field::Perform()
 357 {
 358   if ( myFile.empty() )
 359     return addMessage("File name not set", /*isFatal=*/true ); // 'fatal' means 'bug'
 360   int version = -1, major, minor, release;
 361   if ( MED::GetMEDVersion( myFile, major, minor, release ))
 362     version = major * 10 + minor;
 363
 364   MED::PWrapper medFile = MED::CrWrapperW( myFile, version );
 365   return this->PerformInternal(medFile);
 366 }
 367
 368 /*!
 369  * Writes a field to a chunck of memory
 370  */
 371 Driver_Mesh::Status DriverMED_W_Field_Mem::Perform()
 372 {
 373   Driver_Mesh::Status status = Driver_Mesh::DRS_OK;
 374   bool isClosed(false);
 375   void *ptr(_data->getPointer());
 376   std::size_t sz(_data->getNumberOfTuples());
 377   _data->accessToMemArray().setSpecificDeallocator(nullptr);
 378   _data->useArray(nullptr,false,MEDCoupling::DeallocType::C_DEALLOC,0,1);
 379   {// let braces to flush (call of MED::PWrapper myMed destructor)
 380     MED::TMemFile *tfileInst = new MED::TMemFile(&ptr,&sz,&isClosed);
 381     MED::PWrapper myMed = MED::CrWrapperW(myFile, -1, tfileInst);
 382     status = this->PerformInternal(myMed);
 383   }
 384   if(!isClosed)
 385     EXCEPTION(std::runtime_error, "TFTMemFile destructor : on destruction file has not been closed properly -> chunk of memory data may be invalid !");
 386   _data = MEDCoupling::DataArrayByte::New();
 387   _data->useArray(reinterpret_cast<char *>(ptr),true,MEDCoupling::DeallocType::C_DEALLOC,sz,1);
 388   return status;
 389 }
 390
 391 namespace DriverMED // Implementation of functions declared in DriverMED.hxx
 392 {
 393   //================================================================================
 394   /*!
 395    * Returns a vector containing MED::EGeometrieElement for each SMDSAbs_EntityType
 396    */
 397   //================================================================================
 398
 399   const std::vector< MED::EGeometrieElement >& getMedTypesVec()
 400   {
 401     static std::vector< MED::EGeometrieElement > theVec;
 402     if ( theVec.empty() )
 403     {
 404       theVec.resize( SMDSEntity_Last, MED::eAllGeoType );
 405       theVec[ SMDSEntity_Node               ] = MED::eNONE    ;
 406       theVec[ SMDSEntity_0D                 ] = MED::ePOINT1  ;
 407       theVec[ SMDSEntity_Edge               ] = MED::eSEG2    ;
 408       theVec[ SMDSEntity_Quad_Edge          ] = MED::eSEG3    ;
 409       theVec[ SMDSEntity_Triangle           ] = MED::eTRIA3   ;
 410       theVec[ SMDSEntity_Quad_Triangle      ] = MED::eTRIA6   ;
 411       theVec[ SMDSEntity_BiQuad_Triangle    ] = MED::eTRIA7   ;
 412       theVec[ SMDSEntity_Quadrangle         ] = MED::eQUAD4   ;
 413       theVec[ SMDSEntity_Quad_Quadrangle    ] = MED::eQUAD8   ;
 414       theVec[ SMDSEntity_BiQuad_Quadrangle  ] = MED::eQUAD9   ;
 415       theVec[ SMDSEntity_Polygon            ] = MED::ePOLYGONE;
 416       //theVec[ SMDSEntity_Quad_Polygon       ] = MED::ePOLYGONE; // !!
 417       theVec[ SMDSEntity_Tetra              ] = MED::eTETRA4  ;
 418       theVec[ SMDSEntity_Quad_Tetra         ] = MED::eTETRA10 ;
 419       theVec[ SMDSEntity_Pyramid            ] = MED::ePYRA5   ;
 420       theVec[ SMDSEntity_Quad_Pyramid       ] = MED::ePYRA13  ;
 421       theVec[ SMDSEntity_Hexa               ] = MED::eHEXA8   ;
 422       theVec[ SMDSEntity_Quad_Hexa          ] = MED::eHEXA20  ;
 423       theVec[ SMDSEntity_TriQuad_Hexa       ] = MED::eHEXA27  ;
 424       theVec[ SMDSEntity_Penta              ] = MED::ePENTA6  ;
 425       theVec[ SMDSEntity_Quad_Penta         ] = MED::ePENTA15 ;
 426       theVec[ SMDSEntity_BiQuad_Penta       ] = MED::ePENTA18 ;
 427       theVec[ SMDSEntity_Hexagonal_Prism    ] = MED::eOCTA12  ;
 428       theVec[ SMDSEntity_Polyhedra          ] = MED::ePOLYEDRE;
 429       //theVec[ SMDSEntity_Quad_Polyhedra     ] = MED::ePOLYEDRE; // !!
 430       theVec[ SMDSEntity_Ball               ] = MED::eBALL    ;
 431     }
 432     return theVec;
 433   }
 434
 435   //================================================================================
 436   /*!
 437    * Returns MED element geom type (MED::EGeometrieElement) by SMDS type
 438    */
 439   //================================================================================
 440
 441   int GetMedGeoType( SMDSAbs_EntityType smdsType )
 442   {
 443     return getMedTypesVec()[ smdsType ];
 444   }
 445
 446   //================================================================================
 447   /*!
 448    * Returns SMDS element geom type by MED type (MED::EGeometrieElement)
 449    */
 450   //================================================================================
 451
 452   SMDSAbs_EntityType GetSMDSType( int medType )
 453   {
 454     const std::vector< MED::EGeometrieElement >& theVec = getMedTypesVec();
 455
 456     std::vector< MED::EGeometrieElement >::const_iterator i =
 457       std::find( theVec.begin(), theVec.end(), medType );
 458
 459     return SMDSAbs_EntityType( std::distance( theVec.begin(), i ));
 460   }
 461 }