src/DriverMED/DriverMED_W_Field.cxx

   1 // Copyright (C) 2007-2016  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_Utilities.hxx"
  33 #include "MED_Wrapper.hxx"
  34 #include "SMDS_IteratorOnIterators.hxx"
  35 #include "SMDS_MeshElement.hxx"
  36 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
  37 #include "SMDS_SetIterator.hxx"
  38 #include "SMESHDS_Mesh.hxx"
  39
  40 //================================================================================
  41 /*!
  42  * \brief Constructor
  43  */
  44 //================================================================================
  45
  46 DriverMED_W_Field::DriverMED_W_Field():
  47   //_medFileID( -1 ),
  48   _elemType( SMDSAbs_All ),
  49   _dt( -1 ),
  50   _it( -1 )
  51 {
  52 }
  53
  54 //================================================================================
  55 /*!
  56  * \brief Sets basic data
  57  *  \param [in] mesh - supporting mesh
  58  *  \param [in] fieldName - name of a field
  59  *  \param [in] type - type of supporting elements
  60  *  \param [in] nbComps - number of components
  61  *  \param [in] isIntData - type of data: double or integer
  62  */
  63 //================================================================================
  64
  65 bool DriverMED_W_Field::Set(SMESHDS_Mesh *      mesh,
  66                             const std::string & fieldName,
  67                             SMDSAbs_ElementType type,
  68                             const int           nbComps,
  69                             const bool          isIntData)
  70 {
  71   _fieldName = fieldName;
  72   _compNames.resize( nbComps, "" );
  73
  74   if ( type == SMDSAbs_All )
  75   {
  76     if ( mesh->NbVolumes() > 0 )
  77       type = SMDSAbs_Volume;
  78     else if ( mesh->NbFaces() > 0 )
  79       type = SMDSAbs_Face;
  80     else if ( mesh->NbEdges() > 0 )
  81       type = SMDSAbs_Edge;
  82     else
  83       type = SMDSAbs_Node;
  84   }
  85   if ( myMesh != mesh )
  86   {
  87     _nbElemsByGeom.clear();
  88     for ( int iG = 0; iG < SMDSEntity_Last; ++iG )
  89       _elemsByGeom[iG].clear();
  90     SetMesh( mesh );
  91   }
  92
  93   // find out "MED order" of elements - sort elements by geom type
  94   int nbElems;
  95   if ( _nbElemsByGeom.empty() || _elemType != type )
  96   {
  97     _elemType = type;
  98     _nbElemsByGeom.resize( 1, std::make_pair( SMDSEntity_Last, 0 ));
  99
 100     // count nb of elems of each geometry
 101     for ( int iG = 0; iG < SMDSEntity_Last; ++iG )
 102     {
 103       SMDSAbs_EntityType  geom = (SMDSAbs_EntityType) iG;
 104       SMDSAbs_ElementType t = SMDS_MeshCell::toSmdsType( geom );
 105       if ( t != _elemType ) continue;
 106
 107       nbElems = mesh->GetMeshInfo().NbElements( geom );
 108       if ( nbElems < 1 ) continue;
 109
 110       _nbElemsByGeom.push_back( std::make_pair( geom, nbElems + _nbElemsByGeom.back().second ));
 111     }
 112     // add nodes of missing 0D elements on VERTEXes
 113     if ( _addODOnVertices && _elemType == SMDSAbs_0DElement )
 114     {
 115       std::vector< const SMDS_MeshElement* >& nodes = _elemsByGeom[SMDSEntity_Node];
 116       if ( nodes.empty() )
 117         DriverMED_W_SMESHDS_Mesh::getNodesOfMissing0DOnVert( myMesh, nodes );
 118       if ( !nodes.empty() )
 119       {
 120         if ( _nbElemsByGeom.size() == 1 )
 121           _nbElemsByGeom.push_back( std::make_pair( SMDSEntity_0D, 0));
 122         _nbElemsByGeom.push_back( std::make_pair( SMDSEntity_Node,
 123                                              nodes.size() + _nbElemsByGeom.back().second ));
 124       }
 125     }
 126
 127     // sort elements by their geometry
 128     int iGeoType, nbGeomTypes = _nbElemsByGeom.size() - 1;
 129     if ( nbGeomTypes > 1 )
 130     {
 131       for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 132       {
 133         iGeoType = _nbElemsByGeom[iG].first;
 134         nbElems  = _nbElemsByGeom[iG].second - _nbElemsByGeom[iG-1].second;
 135         _elemsByGeom[ iGeoType ].reserve( nbElems );
 136       }
 137       iGeoType = _nbElemsByGeom[1].first; // for missing 0D
 138       if ( _elemsByGeom[ iGeoType ].empty() )
 139       {
 140         nbElems = mesh->GetMeshInfo().NbElements( _elemType );
 141         SMDS_ElemIteratorPtr eIt = mesh->elementsIterator( _elemType );
 142         for ( int iE = 0; iE < nbElems && eIt->more(); ++iE )
 143         {
 144           const SMDS_MeshElement* e = eIt->next();
 145           _elemsByGeom[ e->GetEntityType() ].push_back( e );
 146         }
 147       }
 148     }
 149   }
 150   _intValues.clear();
 151   _dblValues.clear();
 152
 153   // allocate memory for values
 154   nbElems = _nbElemsByGeom.empty() ? 0 : _nbElemsByGeom.back().second;
 155   if ( isIntData )
 156     _intValues.reserve( nbElems * nbComps );
 157   else
 158     _dblValues.reserve( nbElems * nbComps );
 159
 160   return nbElems * nbComps;
 161 }
 162
 163 //================================================================================
 164 /*!
 165  * \brief Set a name of a component countered from zero
 166  */
 167 //================================================================================
 168
 169 void DriverMED_W_Field::SetCompName(const int iComp, const char* name)
 170 {
 171   if ( (int)_compNames.size() <= iComp )
 172     _compNames.resize( iComp + 1 );
 173   _compNames[ iComp ] = name;
 174 }
 175
 176 //================================================================================
 177 /*!
 178  * \brief Sets numdt and numit field features. Call this fun before AddValue()!
 179  */
 180 //================================================================================
 181
 182 void DriverMED_W_Field::SetDtIt(const int dt, const int it)
 183 {
 184   _dt = dt;
 185   _it = it;
 186   _intValues.clear();
 187   _dblValues.clear();
 188 }
 189
 190 //================================================================================
 191 /*!
 192  * \brief Adds a float field value
 193  */
 194 //================================================================================
 195
 196 void DriverMED_W_Field::AddValue( double val )
 197 {
 198   _dblValues.push_back( val );
 199 }
 200
 201 //================================================================================
 202 /*!
 203  * \brief Adds an integer field value
 204  */
 205 //================================================================================
 206
 207 void DriverMED_W_Field::AddValue( int    val )
 208 {
 209   _intValues.push_back( val );
 210 }
 211
 212 //================================================================================
 213 /*!
 214  * Returns elements in the order they are written in MED file
 215  */
 216 //================================================================================
 217
 218 SMDS_ElemIteratorPtr DriverMED_W_Field::GetOrderedElems()
 219 {
 220   if ( _nbElemsByGeom.size() < 2 )
 221     return SMDS_ElemIteratorPtr();
 222
 223   if ( _nbElemsByGeom.size() == 2 )
 224     // sole geom type of elements
 225     return myMesh->elementsIterator( _elemType );
 226
 227   std::vector< SMDS_ElemIteratorPtr > iterVec( _nbElemsByGeom.size()-1 );
 228   for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 229   {
 230     int iGeoType = _nbElemsByGeom[ iG ].first;
 231     iterVec[ iG-1 ] = SMDS_ElemIteratorPtr
 232       ( new SMDS_ElementVectorIterator( _elemsByGeom[ iGeoType ].begin(),
 233                                         _elemsByGeom[ iGeoType ].end() ));
 234   }
 235   typedef SMDS_IteratorOnIterators
 236     < const SMDS_MeshElement *, std::vector< SMDS_ElemIteratorPtr > > TItIterator;
 237   return SMDS_ElemIteratorPtr( new TItIterator( iterVec ));
 238 }
 239
 240 //================================================================================
 241 /*!
 242  * Writes a field to the file
 243  */
 244 //================================================================================
 245
 246 Driver_Mesh::Status DriverMED_W_Field::Perform()
 247 {
 248   if ( myFile.empty() )
 249     return addMessage("File name not set", /*isFatal=*/true ); // 'fatal' means 'bug'
 250   if ( myMeshId < 0 && myMeshName.empty() )
 251     return addMessage("Mesh in file not specified", /*isFatal=*/true );
 252   if ( _nbElemsByGeom.size() < 2 )
 253     return addMessage("No values to write", /*isFatal=*/false );
 254   if ( !myMesh )
 255     return addMessage("Supporting mesh not set", /*isFatal=*/true );
 256
 257   MED::PWrapper medFile = MED::CrWrapper( myFile, MED::eV2_2 );
 258   MED::PMeshInfo meshInfo;
 259   if ( myMeshId > 0 )
 260   {
 261     meshInfo = medFile->GetPMeshInfo( myMeshId );
 262   }
 263   else
 264   {
 265     // look for a mesh by name
 266     int aNbMeshes = medFile->GetNbMeshes();
 267     for ( int iMesh = aNbMeshes; iMesh > 0 && myMeshId < 1; --iMesh )
 268     {
 269       meshInfo = medFile->GetPMeshInfo( iMesh );
 270       if ( !meshInfo || meshInfo->GetName() != myMeshName )
 271         meshInfo.reset();
 272       else
 273         myMeshId = iMesh;
 274     }
 275   }
 276   if (( !meshInfo ) ||
 277       ( !myMeshName.empty() && meshInfo->GetName() != myMeshName ))
 278   {
 279     myMeshId = -1;
 280     return addMessage("Specified mesh not found in the file", /*isFatal=*/true );
 281   }
 282
 283   // create a field
 284   MED::ETypeChamp  dataType = _dblValues.empty() ? MED::eINT : MED::eFLOAT64;
 285   MED::PFieldInfo fieldInfo = medFile->CrFieldInfo( meshInfo,
 286                                                     _compNames.size(),
 287                                                     dataType );
 288   fieldInfo->SetName( _fieldName );
 289   for ( size_t iC = 0; iC < _compNames.size(); ++iC )
 290   {
 291     fieldInfo->SetCompName( iC, _compNames[ iC ]);
 292     fieldInfo->SetUnitName( iC, "");
 293   }
 294   if ( _compNames.size() > 1 )
 295   {
 296     for ( size_t i = 0; i < fieldInfo->myCompNames.size()-1; ++i )
 297       if ( !fieldInfo->myCompNames[i] )
 298         fieldInfo->myCompNames[i] = ' ';
 299   }
 300   medFile->SetFieldInfo( fieldInfo );
 301
 302   // specific treatment of added 0D elements
 303   if ( _nbElemsByGeom.size()   == 3 &&
 304        _nbElemsByGeom[1].first == SMDSEntity_0D )
 305   {
 306     _nbElemsByGeom[1].second += _nbElemsByGeom[2].second;
 307     _nbElemsByGeom.resize( 2 );
 308   }
 309
 310   // create a time stamp
 311   MED::TGeom2Size type2nb;
 312   for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 313   {
 314     SMDSAbs_EntityType    smdsType = _nbElemsByGeom[iG].first;
 315     MED::EGeometrieElement medType = (MED::EGeometrieElement) DriverMED::GetMedGeoType( smdsType );
 316     int                    nbElems = _nbElemsByGeom[iG].second - _nbElemsByGeom[iG-1].second;
 317     type2nb.insert( std::make_pair( medType, nbElems ));
 318   }
 319
 320   MED::EEntiteMaillage       entity = ( _elemType == SMDSAbs_Node ? MED::eNOEUD : MED::eMAILLE );
 321   MED::PTimeStampInfo timeStampInfo = medFile->CrTimeStampInfo( fieldInfo, entity, type2nb );
 322   timeStampInfo->myNumDt  = _dt;
 323   timeStampInfo->myNumOrd = _it;
 324
 325   MED::PTimeStampValueBase  timeStampVal = medFile->CrTimeStampValue( timeStampInfo, dataType );
 326   MED::PFloatTimeStampValue timeStampFltVal = timeStampVal;
 327   MED::PIntTimeStampValue   timeStampIntVal = timeStampVal;
 328
 329   // set values
 330   int iVal = 0;
 331   MED::TFloat* ptrDbl = 0;
 332   MED::TInt*   ptrInt = 0;
 333   for ( size_t iG = 1; iG < _nbElemsByGeom.size(); ++iG )
 334   {
 335     SMDSAbs_EntityType    smdsType = _nbElemsByGeom[iG].first;
 336     MED::EGeometrieElement medType = (MED::EGeometrieElement) DriverMED::GetMedGeoType( smdsType );
 337     int nbElems = ( _nbElemsByGeom[iG].second - _nbElemsByGeom[iG-1].second ) * _compNames.size();
 338     if ( dataType == MED::eFLOAT64 )
 339     {
 340       ptrDbl = timeStampFltVal->GetMeshValue( medType ).GetPointer();
 341       for ( int i = 0; i < nbElems; ++i, ++iVal )
 342         ptrDbl[ i ] = _dblValues[ iVal ];
 343     }
 344     else
 345     {
 346       ptrInt = timeStampIntVal->GetMeshValue( medType ).GetPointer();
 347       for ( int i = 0; i < nbElems; ++i, ++iVal )
 348         ptrInt[ i ] = _intValues[ iVal ];
 349     }
 350   }
 351
 352   // write
 353   medFile->SetTimeStampValue( timeStampVal );
 354
 355   _dblValues.clear();
 356   _intValues.clear();
 357
 358   return DRS_OK;
 359 }
 360
 361 namespace DriverMED // Implemetation of fuctions declared in DriverMED.hxx
 362 {
 363   //================================================================================
 364   /*!
 365    * Returns a vector containing MED::EGeometrieElement for each SMDSAbs_EntityType
 366    */
 367   //================================================================================
 368
 369   const std::vector< MED::EGeometrieElement >& getMedTypesVec()
 370   {
 371     static std::vector< MED::EGeometrieElement > theVec;
 372     if ( theVec.empty() )
 373     {
 374       theVec.resize( SMDSEntity_Last, MED::eAllGeoType );
 375       theVec[ SMDSEntity_Node               ] = MED::eNONE    ;
 376       theVec[ SMDSEntity_0D                 ] = MED::ePOINT1  ;
 377       theVec[ SMDSEntity_Edge               ] = MED::eSEG2    ;
 378       theVec[ SMDSEntity_Quad_Edge          ] = MED::eSEG3    ;
 379       theVec[ SMDSEntity_Triangle           ] = MED::eTRIA3   ;
 380       theVec[ SMDSEntity_Quad_Triangle      ] = MED::eTRIA6   ;
 381       theVec[ SMDSEntity_BiQuad_Triangle    ] = MED::eTRIA7   ;
 382       theVec[ SMDSEntity_Quadrangle         ] = MED::eQUAD4   ;
 383       theVec[ SMDSEntity_Quad_Quadrangle    ] = MED::eQUAD8   ;
 384       theVec[ SMDSEntity_BiQuad_Quadrangle  ] = MED::eQUAD9   ;
 385       theVec[ SMDSEntity_Polygon            ] = MED::ePOLYGONE;
 386       //theVec[ SMDSEntity_Quad_Polygon       ] = MED::ePOLYGONE; // !!
 387       theVec[ SMDSEntity_Tetra              ] = MED::eTETRA4  ;
 388       theVec[ SMDSEntity_Quad_Tetra         ] = MED::eTETRA10 ;
 389       theVec[ SMDSEntity_Pyramid            ] = MED::ePYRA5   ;
 390       theVec[ SMDSEntity_Quad_Pyramid       ] = MED::ePYRA13  ;
 391       theVec[ SMDSEntity_Hexa               ] = MED::eHEXA8   ;
 392       theVec[ SMDSEntity_Quad_Hexa          ] = MED::eHEXA20  ;
 393       theVec[ SMDSEntity_TriQuad_Hexa       ] = MED::eHEXA27  ;
 394       theVec[ SMDSEntity_Penta              ] = MED::ePENTA6  ;
 395       theVec[ SMDSEntity_Quad_Penta         ] = MED::ePENTA15 ;
 396       theVec[ SMDSEntity_Hexagonal_Prism    ] = MED::eOCTA12  ;
 397       theVec[ SMDSEntity_Polyhedra          ] = MED::ePOLYEDRE;
 398       //theVec[ SMDSEntity_Quad_Polyhedra     ] = MED::ePOLYEDRE; // !!
 399       theVec[ SMDSEntity_Ball               ] = MED::eBALL    ;
 400     }
 401     return theVec;
 402   }
 403
 404   //================================================================================
 405   /*!
 406    * Returns MED element geom type (MED::EGeometrieElement) by SMDS type
 407    */
 408   //================================================================================
 409
 410   int GetMedGeoType( SMDSAbs_EntityType smdsType )
 411   {
 412     return getMedTypesVec()[ smdsType ];
 413   }
 414
 415   //================================================================================
 416   /*!
 417    * Returns SMDS element geom type by MED type (MED::EGeometrieElement)
 418    */
 419   //================================================================================
 420
 421   SMDSAbs_EntityType GetSMDSType( int medType )
 422   {
 423     const std::vector< MED::EGeometrieElement >& theVec = getMedTypesVec();
 424
 425     std::vector< MED::EGeometrieElement >::const_iterator i =
 426       std::find( theVec.begin(), theVec.end(), medType );
 427
 428     return SMDSAbs_EntityType( std::distance( theVec.begin(), i ));
 429   }
 430 }