Salome HOME
Update copyright information
[modules/smesh.git] / src / DriverMED / DriverMED_W_SMESHDS_Mesh.cxx
1 // Copyright (C) 2007-2012  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.
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 //  SMESH DriverMED : driver to read and write 'med' files
24 //  File   : DriverMED_W_SMESHDS_Mesh.cxx
25 //  Module : SMESH
26 //
27 #include <sstream>
28
29 #include "DriverMED_W_SMESHDS_Mesh.h"
30 #include "DriverMED_Family.h"
31
32 #include "SMESHDS_Mesh.hxx"
33 #include "SMDS_MeshElement.hxx"
34 #include "SMDS_MeshNode.hxx"
35 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
36
37 #include "utilities.h"
38
39 #include "MED_Utilities.hxx"
40
41 #define _EDF_NODE_IDS_
42 //#define _ELEMENTS_BY_DIM_
43
44 using namespace std;
45 using namespace MED;
46
47
48 DriverMED_W_SMESHDS_Mesh::DriverMED_W_SMESHDS_Mesh():
49   myAllSubMeshes (false),
50   myDoGroupOfNodes (false),
51   myDoGroupOfEdges (false),
52   myDoGroupOfFaces (false),
53   myDoGroupOfVolumes (false)
54 {}
55
56 void DriverMED_W_SMESHDS_Mesh::SetFile(const std::string& theFileName, 
57                                        MED::EVersion theId)
58 {
59   myMed = CrWrapper(theFileName,theId);
60   Driver_SMESHDS_Mesh::SetFile(theFileName);
61 }
62
63 void DriverMED_W_SMESHDS_Mesh::SetFile(const std::string& theFileName)
64 {
65   return SetFile(theFileName,MED::eV2_2);
66 }
67
68 string DriverMED_W_SMESHDS_Mesh::GetVersionString(const MED::EVersion theVersion, int theNbDigits)
69 {
70   TInt majeur, mineur, release;
71   majeur =  mineur = release = 0;
72 //   if ( theVersion == eV2_1 )
73 //     MED::GetVersionRelease<eV2_1>(majeur, mineur, release);
74 //   else
75     MED::GetVersionRelease<eV2_2>(majeur, mineur, release);
76   ostringstream name;
77   if ( theNbDigits > 0 )
78     name << majeur;
79   if ( theNbDigits > 1 )
80     name << "." << mineur;
81   if ( theNbDigits > 2 )
82     name << "." << release;
83   return name.str();
84 }
85
86 void DriverMED_W_SMESHDS_Mesh::AddGroup(SMESHDS_GroupBase* theGroup)
87 {
88   myGroups.push_back(theGroup);
89 }
90
91 void DriverMED_W_SMESHDS_Mesh::AddAllSubMeshes()
92 {
93   myAllSubMeshes = true;
94 }
95
96 void DriverMED_W_SMESHDS_Mesh::AddSubMesh(SMESHDS_SubMesh* theSubMesh, int theID)
97 {
98   mySubMeshes[theID] = theSubMesh;
99 }
100
101 void DriverMED_W_SMESHDS_Mesh::AddGroupOfNodes()
102 {
103   myDoGroupOfNodes = true;
104 }
105
106 void DriverMED_W_SMESHDS_Mesh::AddGroupOfEdges()
107 {
108   myDoGroupOfEdges = true;
109 }
110
111 void DriverMED_W_SMESHDS_Mesh::AddGroupOfFaces()
112 {
113   myDoGroupOfFaces = true;
114 }
115
116 void DriverMED_W_SMESHDS_Mesh::AddGroupOfVolumes()
117 {
118   myDoGroupOfVolumes = true;
119 }
120
121 namespace{
122   typedef double (SMDS_MeshNode::* TGetCoord)() const;
123   typedef const char* TName;
124   typedef const char* TUnit;
125
126   // name length in a mesh must be equal to 16 :
127   //         1234567890123456
128   TName M = "m               ";
129   TName X = "x               ";
130   TName Y = "y               ";
131   TName Z = "z               ";
132
133   TUnit aUnit[3] = {M,M,M};
134
135   // 3 dim
136   TGetCoord aXYZGetCoord[3] = {
137     &SMDS_MeshNode::X, 
138     &SMDS_MeshNode::Y, 
139     &SMDS_MeshNode::Z
140   };
141   TName aXYZName[3] = {X,Y,Z};
142   
143   // 2 dim
144   TGetCoord aXYGetCoord[2] = {
145     &SMDS_MeshNode::X, 
146     &SMDS_MeshNode::Y
147   };
148   TName aXYName[2] = {X,Y};
149
150   TGetCoord aYZGetCoord[2] = {
151     &SMDS_MeshNode::Y, 
152     &SMDS_MeshNode::Z
153   };
154   TName aYZName[2] = {Y,Z};
155
156   TGetCoord aXZGetCoord[2] = {
157     &SMDS_MeshNode::X, 
158     &SMDS_MeshNode::Z
159   };
160   TName aXZName[2] = {X,Z};
161
162   // 1 dim
163   TGetCoord aXGetCoord[1] = {
164     &SMDS_MeshNode::X
165   };
166   TName aXName[1] = {X};
167
168   TGetCoord aYGetCoord[1] = {
169     &SMDS_MeshNode::Y
170   };
171   TName aYName[1] = {Y};
172
173   TGetCoord aZGetCoord[1] = {
174     &SMDS_MeshNode::Z
175   };
176   TName aZName[1] = {Z};
177
178
179   class TCoordHelper{
180     SMDS_NodeIteratorPtr myNodeIter;
181     const SMDS_MeshNode* myCurrentNode;
182     TGetCoord* myGetCoord;
183     TName* myName;
184     TUnit* myUnit;
185   public:
186     TCoordHelper(const SMDS_NodeIteratorPtr& theNodeIter,
187                  TGetCoord* theGetCoord,
188                  TName* theName,
189                  TUnit* theUnit = aUnit):
190       myNodeIter(theNodeIter),
191       myGetCoord(theGetCoord),
192       myName(theName),
193       myUnit(theUnit)
194     {}
195     virtual ~TCoordHelper(){}
196     bool Next(){ 
197       return myNodeIter->more() && 
198         (myCurrentNode = myNodeIter->next());
199     }
200     const SMDS_MeshNode* GetNode(){
201       return myCurrentNode;
202     }
203     MED::TIntVector::value_type GetID(){
204       return myCurrentNode->GetID();
205     }
206     MED::TFloatVector::value_type GetCoord(TInt theCoodId){
207       return (myCurrentNode->*myGetCoord[theCoodId])();
208     }
209     MED::TStringVector::value_type GetName(TInt theDimId){
210       return myName[theDimId];
211     }
212     MED::TStringVector::value_type GetUnit(TInt theDimId){
213       return myUnit[theDimId];
214     }
215   };
216   typedef boost::shared_ptr<TCoordHelper> TCoordHelperPtr;
217
218   //-------------------------------------------------------
219   /*!
220    * \brief Structure describing element type
221    */
222   //-------------------------------------------------------
223   struct TElemTypeData
224   {
225     EEntiteMaillage     _entity;
226     EGeometrieElement   _geomType;
227     TInt                _nbElems;
228     SMDSAbs_ElementType _smdsType;
229
230     TElemTypeData (EEntiteMaillage entity, EGeometrieElement geom, TInt nb, SMDSAbs_ElementType type)
231       : _entity(entity), _geomType(geom), _nbElems( nb ), _smdsType( type ) {}
232   };
233
234
235   typedef NCollection_DataMap< Standard_Address, int > TElemFamilyMap;
236
237   //================================================================================
238   /*!
239    * \brief Fills element to famaly ID map for element type.
240    * Removes all families of anElemType
241    */
242   //================================================================================
243
244   void fillElemFamilyMap( TElemFamilyMap &            anElemFamMap,
245                           list<DriverMED_FamilyPtr> & aFamilies,
246                           const SMDSAbs_ElementType   anElemType)
247   {
248     anElemFamMap.Clear();
249     //anElemFamMap.clear();
250     list<DriverMED_FamilyPtr>::iterator aFamsIter = aFamilies.begin();
251     while ( aFamsIter != aFamilies.end() )
252     {
253       if ((*aFamsIter)->GetType() != anElemType) {
254         aFamsIter++;
255       }
256       else {
257         int aFamId = (*aFamsIter)->GetId();
258         const set<const SMDS_MeshElement *>& anElems = (*aFamsIter)->GetElements();
259         set<const SMDS_MeshElement *>::const_iterator anElemsIter = anElems.begin();
260         for (; anElemsIter != anElems.end(); anElemsIter++)
261         {
262           anElemFamMap.Bind( (Standard_Address)*anElemsIter, aFamId );
263           //anElemFamMap[*anElemsIter] = aFamId;
264         }
265         // remove a family from the list
266         aFamilies.erase( aFamsIter++ );
267       }
268     }
269   }
270
271   //================================================================================
272   /*!
273    * \brief For an element, return family ID found in the map or a default one
274    */
275   //================================================================================
276
277   int getFamilyId( const TElemFamilyMap &  anElemFamMap,
278                    const SMDS_MeshElement* anElement,
279                    const int               aDefaultFamilyId)
280   {
281     if ( anElemFamMap.IsBound( (Standard_Address) anElement ))
282       return anElemFamMap( (Standard_Address) anElement );
283 //     TElemFamilyMap::iterator elem_famNum = anElemFamMap.find( anElement );
284 //     if ( elem_famNum != anElemFamMap.end() )
285 //       return elem_famNum->second;
286     return aDefaultFamilyId;
287   }
288 }
289
290 Driver_Mesh::Status DriverMED_W_SMESHDS_Mesh::Perform()
291 {
292   Status aResult = DRS_OK;
293   if (myMesh->hasConstructionEdges() || myMesh->hasConstructionFaces()) {
294     INFOS("SMDS_MESH with hasConstructionEdges() or hasConstructionFaces() do not supports!!!");
295     return DRS_FAIL;
296   }
297   try {
298     MESSAGE("Perform - myFile : "<<myFile);
299
300     // Creating the MED mesh for corresponding SMDS structure
301     //-------------------------------------------------------
302     string aMeshName;
303     if (myMeshId != -1) {
304       ostringstream aMeshNameStr;
305       aMeshNameStr<<myMeshId;
306       aMeshName = aMeshNameStr.str();
307     } else {
308       aMeshName = myMeshName;
309     }
310
311     // Mesh dimension definition
312     TInt aSpaceDimension;
313     TCoordHelperPtr aCoordHelperPtr;
314     {  
315       bool anIsXDimension = false;
316       bool anIsYDimension = false;
317       bool anIsZDimension = false;
318       {
319         SMDS_NodeIteratorPtr aNodesIter = myMesh->nodesIterator();
320         double aBounds[6];
321         if(aNodesIter->more()){
322           const SMDS_MeshNode* aNode = aNodesIter->next();
323           aBounds[0] = aBounds[1] = aNode->X();
324           aBounds[2] = aBounds[3] = aNode->Y();
325           aBounds[4] = aBounds[5] = aNode->Z();
326         }
327         while(aNodesIter->more()){
328           const SMDS_MeshNode* aNode = aNodesIter->next();
329           aBounds[0] = min(aBounds[0],aNode->X());
330           aBounds[1] = max(aBounds[1],aNode->X());
331           
332           aBounds[2] = min(aBounds[2],aNode->Y());
333           aBounds[3] = max(aBounds[3],aNode->Y());
334           
335           aBounds[4] = min(aBounds[4],aNode->Z());
336           aBounds[5] = max(aBounds[5],aNode->Z());
337         }
338
339         double EPS = 1.0E-7;
340         anIsXDimension = (aBounds[1] - aBounds[0]) + abs(aBounds[1]) + abs(aBounds[0]) > EPS;
341         anIsYDimension = (aBounds[3] - aBounds[2]) + abs(aBounds[3]) + abs(aBounds[2]) > EPS;
342         anIsZDimension = (aBounds[5] - aBounds[4]) + abs(aBounds[5]) + abs(aBounds[4]) > EPS;
343         aSpaceDimension = anIsXDimension + anIsYDimension + anIsZDimension;
344         if(!aSpaceDimension)
345           aSpaceDimension = 3;
346         // PAL16857(SMESH not conform to the MED convention):
347         if ( aSpaceDimension == 2 && anIsZDimension ) // 2D only if mesh is in XOY plane
348           aSpaceDimension = 3;
349         // PAL18941(a saved study with a mesh belong Z is opened and the mesh is belong X)
350         if ( aSpaceDimension == 1 && !anIsXDimension ) {// 1D only if mesh is along OX
351           if ( anIsYDimension ) {
352             aSpaceDimension = 2;
353             anIsXDimension = true;
354           } else {
355             aSpaceDimension = 3;
356           }
357         }
358       }
359
360       SMDS_NodeIteratorPtr aNodesIter = myMesh->nodesIterator(/*idInceasingOrder=*/true);
361       switch(aSpaceDimension){
362       case 3:
363         aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXYZGetCoord,aXYZName));
364         break;
365       case 2:
366         if(anIsXDimension && anIsYDimension)
367           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXYGetCoord,aXYName));
368         if(anIsYDimension && anIsZDimension)
369           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aYZGetCoord,aYZName));
370         if(anIsXDimension && anIsZDimension)
371           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXZGetCoord,aXZName));
372         break;
373       case 1:
374         if(anIsXDimension)
375           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXGetCoord,aXName));
376         if(anIsYDimension)
377           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aYGetCoord,aYName));
378         if(anIsZDimension)
379           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aZGetCoord,aZName));
380         break;
381       }
382     }
383     TInt aMeshDimension = 0;
384     if ( myMesh->NbEdges() > 0 )
385       aMeshDimension = 1;
386     if ( myMesh->NbFaces() > 0 )
387       aMeshDimension = 2;
388     if ( myMesh->NbVolumes() > 0 )
389       aMeshDimension = 3;
390     
391     PMeshInfo aMeshInfo = myMed->CrMeshInfo(aMeshDimension,aSpaceDimension,aMeshName);
392     MESSAGE("Add - aMeshName : "<<aMeshName<<"; "<<aMeshInfo->GetName());
393     myMed->SetMeshInfo(aMeshInfo);
394
395     // Storing SMDS groups and sub-meshes as med families
396     //----------------------------------------------------
397     int myNodesDefaultFamilyId   = 0;
398     int my0DElementsDefaultFamilyId = 0;
399     int myEdgesDefaultFamilyId   = 0;
400     int myFacesDefaultFamilyId   = 0;
401     int myVolumesDefaultFamilyId = 0;
402     int nbNodes   = myMesh->NbNodes();
403     //int nb0DElements = myMesh->Nb0DElements();
404     int nbEdges   = myMesh->NbEdges();
405     int nbFaces   = myMesh->NbFaces();
406     int nbVolumes = myMesh->NbVolumes();
407     if (myDoGroupOfNodes && nbNodes)
408       myNodesDefaultFamilyId = REST_NODES_FAMILY;
409     if (myDoGroupOfEdges && nbEdges)
410       myEdgesDefaultFamilyId = REST_EDGES_FAMILY;
411     if (myDoGroupOfFaces && nbFaces)
412       myFacesDefaultFamilyId = REST_FACES_FAMILY;
413     if (myDoGroupOfVolumes && nbVolumes)
414       myVolumesDefaultFamilyId = REST_VOLUMES_FAMILY;
415
416     MESSAGE("Perform - aFamilyInfo");
417     //cout << " DriverMED_Family::MakeFamilies() " << endl;
418     list<DriverMED_FamilyPtr> aFamilies;
419     if (myAllSubMeshes) {
420       aFamilies = DriverMED_Family::MakeFamilies
421         (myMesh->SubMeshes(), myGroups,
422          myDoGroupOfNodes   && nbNodes,
423          myDoGroupOfEdges   && nbEdges,
424          myDoGroupOfFaces   && nbFaces,
425          myDoGroupOfVolumes && nbVolumes);
426     } else {
427       aFamilies = DriverMED_Family::MakeFamilies
428         (mySubMeshes, myGroups,
429          myDoGroupOfNodes   && nbNodes,
430          myDoGroupOfEdges   && nbEdges,
431          myDoGroupOfFaces   && nbFaces,
432          myDoGroupOfVolumes && nbVolumes);
433     }
434     //cout << " myMed->SetFamilyInfo() " << endl;
435     list<DriverMED_FamilyPtr>::iterator aFamsIter;
436     for (aFamsIter = aFamilies.begin(); aFamsIter != aFamilies.end(); aFamsIter++)
437     {
438       PFamilyInfo aFamilyInfo = (*aFamsIter)->GetFamilyInfo(myMed,aMeshInfo);
439       myMed->SetFamilyInfo(aFamilyInfo);
440     }
441
442     // Storing SMDS nodes to the MED file for the MED mesh
443     //----------------------------------------------------
444 #ifdef _EDF_NODE_IDS_
445     typedef map<TInt,TInt> TNodeIdMap;
446     TNodeIdMap aNodeIdMap;
447 #endif
448     const EModeSwitch   theMode        = eFULL_INTERLACE;
449     const ERepere       theSystem      = eCART;
450     const EBooleen      theIsElemNum   = eVRAI;
451     const EBooleen      theIsElemNames = eFAUX;
452     const EConnectivite theConnMode    = eNOD;
453
454     TInt aNbNodes = myMesh->NbNodes();
455     //cout << " myMed->CrNodeInfo() aNbNodes = " << aNbNodes << endl;
456     PNodeInfo aNodeInfo = myMed->CrNodeInfo(aMeshInfo, aNbNodes,
457                                             theMode, theSystem, theIsElemNum, theIsElemNames);
458
459     //cout << " fillElemFamilyMap( SMDSAbs_Node )" << endl;
460     // find family numbers for nodes
461     TElemFamilyMap anElemFamMap;
462     fillElemFamilyMap( anElemFamMap, aFamilies, SMDSAbs_Node );
463
464     for (TInt iNode = 0; aCoordHelperPtr->Next(); iNode++)
465     {
466       // coordinates
467       TCoordSlice aTCoordSlice = aNodeInfo->GetCoordSlice( iNode );
468       for(TInt iCoord = 0; iCoord < aSpaceDimension; iCoord++){
469         aTCoordSlice[iCoord] = aCoordHelperPtr->GetCoord(iCoord);
470       }
471       // node number
472       int aNodeID = aCoordHelperPtr->GetID();
473       aNodeInfo->SetElemNum( iNode, aNodeID );
474 #ifdef _EDF_NODE_IDS_
475       aNodeIdMap[aNodeID] = iNode+1;
476 #endif
477       // family number
478       const SMDS_MeshNode* aNode = aCoordHelperPtr->GetNode();
479       int famNum = getFamilyId( anElemFamMap, aNode, myNodesDefaultFamilyId );
480       aNodeInfo->SetFamNum( iNode, famNum );
481     }
482     anElemFamMap.Clear();
483
484     // coordinate names and units
485     for (TInt iCoord = 0; iCoord < aSpaceDimension; iCoord++) {
486       aNodeInfo->SetCoordName( iCoord, aCoordHelperPtr->GetName(iCoord));
487       aNodeInfo->SetCoordUnit( iCoord, aCoordHelperPtr->GetUnit(iCoord));
488     }
489
490     //cout << " SetNodeInfo(aNodeInfo)" << endl;
491     MESSAGE("Perform - aNodeInfo->GetNbElem() = "<<aNbNodes);
492     myMed->SetNodeInfo(aNodeInfo);
493     aNodeInfo.reset(); // free memory used for arrays
494
495
496     // Storing SMDS elements to the MED file for the MED mesh
497     //-------------------------------------------------------
498     // Write one element type at once in order to minimize memory usage (PAL19276)
499
500     const SMDS_MeshInfo& nbElemInfo = myMesh->GetMeshInfo();
501
502     // poly elements are not supported by med-2.1
503     bool polyTypesSupported = myMed->CrPolygoneInfo(aMeshInfo,eMAILLE,ePOLYGONE,0,0);
504     TInt nbPolygonNodes = 0, nbPolyhedronNodes = 0, nbPolyhedronFaces = 0;
505
506     // collect info on all geom types
507
508     list< TElemTypeData > aTElemTypeDatas;
509
510     EEntiteMaillage anEntity = eMAILLE;
511 #ifdef _ELEMENTS_BY_DIM_
512     anEntity = eNOEUD_ELEMENT;
513 #endif
514     aTElemTypeDatas.push_back(TElemTypeData(anEntity,
515                                             ePOINT1,
516                                             nbElemInfo.Nb0DElements(),
517                                             SMDSAbs_0DElement));
518 #ifdef _ELEMENTS_BY_DIM_
519     anEntity = eARETE;
520 #endif
521     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
522                                             eSEG2,
523                                             nbElemInfo.NbEdges( ORDER_LINEAR ),
524                                             SMDSAbs_Edge));
525     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
526                                             eSEG3,
527                                             nbElemInfo.NbEdges( ORDER_QUADRATIC ),
528                                             SMDSAbs_Edge));
529 #ifdef _ELEMENTS_BY_DIM_
530     anEntity = eFACE;
531 #endif
532     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
533                                             eTRIA3,
534                                             nbElemInfo.NbTriangles( ORDER_LINEAR ),
535                                             SMDSAbs_Face));
536     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
537                                             eTRIA6,
538                                             nbElemInfo.NbTriangles( ORDER_QUADRATIC ),
539                                             SMDSAbs_Face));
540     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
541                                             eQUAD4,
542                                             nbElemInfo.NbQuadrangles( ORDER_LINEAR ),
543                                             SMDSAbs_Face));
544     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
545                                              eQUAD8,
546                                              nbElemInfo.NbQuadrangles( ORDER_QUADRATIC ) -
547                                              nbElemInfo.NbBiQuadQuadrangles(),
548                                              SMDSAbs_Face));
549     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
550                                              eQUAD9,
551                                              nbElemInfo.NbBiQuadQuadrangles(),
552                                              SMDSAbs_Face));
553     if ( polyTypesSupported ) {
554       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
555                                                ePOLYGONE,
556                                                nbElemInfo.NbPolygons(),
557                                                SMDSAbs_Face));
558       // we need one more loop on poly elements to count nb of their nodes
559       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
560                                                ePOLYGONE,
561                                                nbElemInfo.NbPolygons(),
562                                                SMDSAbs_Face));
563     }
564 #ifdef _ELEMENTS_BY_DIM_
565     anEntity = eMAILLE;
566 #endif
567     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
568                                             eTETRA4,
569                                             nbElemInfo.NbTetras( ORDER_LINEAR ),
570                                             SMDSAbs_Volume));
571     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
572                                             eTETRA10,
573                                             nbElemInfo.NbTetras( ORDER_QUADRATIC ),
574                                             SMDSAbs_Volume));
575     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
576                                             ePYRA5,
577                                             nbElemInfo.NbPyramids( ORDER_LINEAR ),
578                                             SMDSAbs_Volume));
579     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
580                                             ePYRA13,
581                                             nbElemInfo.NbPyramids( ORDER_QUADRATIC ),
582                                             SMDSAbs_Volume));
583     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
584                                             ePENTA6,
585                                             nbElemInfo.NbPrisms( ORDER_LINEAR ),
586                                             SMDSAbs_Volume));
587     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
588                                              ePENTA15,
589                                              nbElemInfo.NbPrisms( ORDER_QUADRATIC ),
590                                              SMDSAbs_Volume));
591     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
592                                              eHEXA8,
593                                              nbElemInfo.NbHexas( ORDER_LINEAR ),
594                                              SMDSAbs_Volume));
595     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
596                                              eHEXA20,
597                                              nbElemInfo.NbHexas( ORDER_QUADRATIC )-
598                                              nbElemInfo.NbTriQuadHexas(),
599                                              SMDSAbs_Volume));
600     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
601                                              eHEXA27,
602                                              nbElemInfo.NbTriQuadHexas(),
603                                              SMDSAbs_Volume));
604     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
605                                              eOCTA12,
606                                              nbElemInfo.NbHexPrisms(),
607                                              SMDSAbs_Volume));
608     if ( polyTypesSupported ) {
609       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
610                                                ePOLYEDRE,
611                                                nbElemInfo.NbPolyhedrons(),
612                                                SMDSAbs_Volume));
613       // we need one more loop on poly elements to count nb of their nodes
614       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
615                                                ePOLYEDRE,
616                                                nbElemInfo.NbPolyhedrons(),
617                                                SMDSAbs_Volume));
618     }
619
620     vector< bool > isElemFamMapBuilt( SMDSAbs_NbElementTypes, false );
621
622     // loop on all geom types of elements
623
624     list< TElemTypeData >::iterator aElemTypeData = aTElemTypeDatas.begin();
625     for ( ; aElemTypeData != aTElemTypeDatas.end(); ++aElemTypeData )
626     {
627       if ( aElemTypeData->_nbElems == 0 )
628         continue;
629
630       int defaultFamilyId = 0;
631       switch ( aElemTypeData->_smdsType ) {
632       case SMDSAbs_0DElement:
633         defaultFamilyId = my0DElementsDefaultFamilyId;
634         break;
635       case SMDSAbs_Edge:
636         defaultFamilyId = myEdgesDefaultFamilyId;
637         break;
638       case SMDSAbs_Face:
639         defaultFamilyId = myFacesDefaultFamilyId;
640         break;
641       case SMDSAbs_Volume:
642         defaultFamilyId = myVolumesDefaultFamilyId;
643         break;
644       default:
645         continue;
646       }
647
648       // iterator on elements of a current type
649       SMDS_ElemIteratorPtr elemIterator = myMesh->elementsIterator( aElemTypeData->_smdsType );
650       if ( !elemIterator->more()) continue;
651       int iElem = 0;
652
653       // Treat POLYGONs
654       // ---------------
655       if ( aElemTypeData->_geomType == ePOLYGONE )
656       {
657         if ( nbPolygonNodes == 0 ) {
658           // Count nb of nodes
659           while ( const SMDS_MeshElement* anElem = elemIterator->next() ) {
660             if ( anElem->IsPoly() ) {
661               nbPolygonNodes += anElem->NbNodes();
662               if ( ++iElem == aElemTypeData->_nbElems )
663                 break;
664             }
665           }
666         }
667         else {
668           // Store in med file
669           PPolygoneInfo aPolygoneInfo = myMed->CrPolygoneInfo(aMeshInfo,
670                                                               aElemTypeData->_entity,
671                                                               aElemTypeData->_geomType,
672                                                               aElemTypeData->_nbElems,
673                                                               nbPolygonNodes,
674                                                               theConnMode, theIsElemNum,
675                                                               theIsElemNames);
676           TElemNum & index = *(aPolygoneInfo->myIndex.get());
677           index[0] = 1;
678
679           while ( const SMDS_MeshElement* anElem = elemIterator->next() )
680           {
681             if ( !anElem->IsPoly() )
682               continue;
683
684             // index
685             TInt aNbNodes = anElem->NbNodes();
686             index[ iElem+1 ] = index[ iElem ] + aNbNodes;
687
688             // connectivity
689             TConnSlice aTConnSlice = aPolygoneInfo->GetConnSlice( iElem );
690             for(TInt iNode = 0; iNode < aNbNodes; iNode++) {
691               const SMDS_MeshElement* aNode = anElem->GetNode( iNode );
692 #ifdef _EDF_NODE_IDS_
693               aTConnSlice[ iNode ] = aNodeIdMap[aNode->GetID()];
694 #else
695               aTConnSlice[ iNode ] = aNode->GetID();
696 #endif
697             }
698             // element number
699             aPolygoneInfo->SetElemNum( iElem, anElem->GetID() );
700
701             // family number
702             int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
703             aPolygoneInfo->SetFamNum( iElem, famNum );
704
705             if ( ++iElem == aPolygoneInfo->GetNbElem() )
706               break;
707           }
708           myMed->SetPolygoneInfo(aPolygoneInfo);
709         }
710
711       } 
712
713       // Treat POLYEDREs
714       // ----------------
715       else if (aElemTypeData->_geomType == ePOLYEDRE )
716       {
717         if ( nbPolyhedronNodes == 0 ) {
718           // Count nb of nodes
719           while ( const SMDS_MeshElement* anElem = elemIterator->next() ) {
720               const SMDS_VtkVolume *aPolyedre = dynamic_cast<const SMDS_VtkVolume*>(anElem);
721               if ( aPolyedre && aPolyedre->IsPoly()) {
722               nbPolyhedronNodes += aPolyedre->NbNodes();
723               nbPolyhedronFaces += aPolyedre->NbFaces();
724               if ( ++iElem == aElemTypeData->_nbElems )
725                 break;
726             }
727           }
728         }
729         else {
730           // Store in med file
731           PPolyedreInfo aPolyhInfo = myMed->CrPolyedreInfo(aMeshInfo,
732                                                            aElemTypeData->_entity,
733                                                            aElemTypeData->_geomType,
734                                                            aElemTypeData->_nbElems,
735                                                            nbPolyhedronFaces+1,
736                                                            nbPolyhedronNodes,
737                                                            theConnMode,
738                                                            theIsElemNum,
739                                                            theIsElemNames);
740           TElemNum & index = *(aPolyhInfo->myIndex.get());
741           TElemNum & faces = *(aPolyhInfo->myFaces.get());
742           TElemNum & conn  = *(aPolyhInfo->myConn.get());
743           index[0] = 1;
744           faces[0] = 1;
745
746           TInt iFace = 0, iNode = 0;
747           while ( const SMDS_MeshElement* anElem = elemIterator->next() )
748           {
749             const SMDS_VtkVolume *aPolyedre = dynamic_cast<const SMDS_VtkVolume*>(anElem);
750             if ( !aPolyedre )
751               continue;
752             if ( !aPolyedre->IsPoly() )
753               continue;
754             // index
755             TInt aNbFaces = aPolyedre->NbFaces();
756             index[ iElem+1 ] = index[ iElem ] + aNbFaces;
757
758             // face index
759             for (TInt f = 1; f <= aNbFaces; ++f, ++iFace ) {
760               int aNbFaceNodes = aPolyedre->NbFaceNodes( f );
761               faces[ iFace+1 ] = faces[ iFace ] + aNbFaceNodes;
762             }
763             // connectivity
764             SMDS_ElemIteratorPtr nodeIt = anElem->nodesIterator();
765             while ( nodeIt->more() ) {
766               const SMDS_MeshElement* aNode = nodeIt->next();
767 #ifdef _EDF_NODE_IDS_
768               conn[ iNode ] = aNodeIdMap[aNode->GetID()];
769 #else
770               conn[ iNode ] = aNode->GetID();
771 #endif
772               ++iNode;
773             }
774             // element number
775             aPolyhInfo->SetElemNum( iElem, anElem->GetID() );
776
777             // family number
778             int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
779             aPolyhInfo->SetFamNum( iElem, famNum );
780
781             if ( ++iElem == aPolyhInfo->GetNbElem() )
782               break;
783           }
784           myMed->SetPolyedreInfo(aPolyhInfo);
785         }
786       } // if (aElemTypeData->_geomType == ePOLYEDRE )
787
788       else
789       {
790         // Treat standard types
791         // ---------------------
792
793         // allocate data arrays
794         PCellInfo aCellInfo = myMed->CrCellInfo( aMeshInfo,
795                                                  aElemTypeData->_entity,
796                                                  aElemTypeData->_geomType,
797                                                  aElemTypeData->_nbElems,
798                                                  theConnMode,
799                                                  theIsElemNum,
800                                                  theIsElemNames);
801         // build map of family numbers for this type
802         if ( !isElemFamMapBuilt[ aElemTypeData->_smdsType ])
803         {
804           //cout << " fillElemFamilyMap()" << endl;
805           fillElemFamilyMap( anElemFamMap, aFamilies, aElemTypeData->_smdsType );
806           isElemFamMapBuilt[ aElemTypeData->_smdsType ] = true;
807         }
808
809         TInt aNbNodes = MED::GetNbNodes(aElemTypeData->_geomType);
810         while ( const SMDS_MeshElement* anElem = elemIterator->next() )
811         {
812           if ( anElem->NbNodes() != aNbNodes || anElem->IsPoly() )
813             continue; // other geometry
814
815           // connectivity
816           TConnSlice aTConnSlice = aCellInfo->GetConnSlice( iElem );
817           for (TInt iNode = 0; iNode < aNbNodes; iNode++) {
818             const SMDS_MeshElement* aNode = anElem->GetNode( iNode );
819 #ifdef _EDF_NODE_IDS_
820             aTConnSlice[ iNode ] = aNodeIdMap[aNode->GetID()];
821 #else
822             aTConnSlice[ iNode ] = aNode->GetID();
823 #endif
824           }
825           // element number
826           aCellInfo->SetElemNum( iElem, anElem->GetID() );
827
828           // family number
829           int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
830           aCellInfo->SetFamNum( iElem, famNum );
831
832           if ( ++iElem == aCellInfo->GetNbElem() )
833             break;
834         }
835         // store data in a file
836         myMed->SetCellInfo(aCellInfo);
837       }
838
839     } // loop on geom types
840
841
842   }
843   catch(const std::exception& exc) {
844     INFOS("The following exception was caught:\n\t"<<exc.what());
845     throw;
846   }
847   catch(...) {
848     INFOS("Unknown exception was caught !!!");
849     throw;
850   }
851
852   myMeshId = -1;
853   myGroups.clear();
854   mySubMeshes.clear();
855   return aResult;
856 }