Salome HOME
23514: EDF 16031 - SMESH freezes
[modules/smesh.git] / 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_BiQuad_Penta       ] = MED::ePENTA18 ;
397       theVec[ SMDSEntity_Hexagonal_Prism    ] = MED::eOCTA12  ;
398       theVec[ SMDSEntity_Polyhedra          ] = MED::ePOLYEDRE;
399       //theVec[ SMDSEntity_Quad_Polyhedra     ] = MED::ePOLYEDRE; // !!
400       theVec[ SMDSEntity_Ball               ] = MED::eBALL    ;
401     }
402     return theVec;
403   }
404
405   //================================================================================
406   /*!
407    * Returns MED element geom type (MED::EGeometrieElement) by SMDS type
408    */
409   //================================================================================
410
411   int GetMedGeoType( SMDSAbs_EntityType smdsType )
412   {
413     return getMedTypesVec()[ smdsType ];
414   }
415
416   //================================================================================
417   /*!
418    * Returns SMDS element geom type by MED type (MED::EGeometrieElement)
419    */
420   //================================================================================
421
422   SMDSAbs_EntityType GetSMDSType( int medType )
423   {
424     const std::vector< MED::EGeometrieElement >& theVec = getMedTypesVec();
425
426     std::vector< MED::EGeometrieElement >::const_iterator i =
427       std::find( theVec.begin(), theVec.end(), medType );
428
429     return SMDSAbs_EntityType( std::distance( theVec.begin(), i ));
430   }
431 }