src/DriverMED/DriverMED_W_Field.cxx

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