Salome HOME
Merge from V6_main_20120808 08Aug12
[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 myBallsDefaultFamilyId      = 0;
400     int myEdgesDefaultFamilyId      = 0;
401     int myFacesDefaultFamilyId      = 0;
402     int myVolumesDefaultFamilyId    = 0;
403     int nbNodes      = myMesh->NbNodes();
404     int nb0DElements = myMesh->Nb0DElements();
405     int nbBalls      = myMesh->NbBalls();
406     int nbEdges      = myMesh->NbEdges();
407     int nbFaces      = myMesh->NbFaces();
408     int nbVolumes    = myMesh->NbVolumes();
409     if (myDoGroupOfNodes && nbNodes) myNodesDefaultFamilyId = REST_NODES_FAMILY;
410     if (myDoGroupOfEdges && nbEdges) myEdgesDefaultFamilyId = REST_EDGES_FAMILY;
411     if (myDoGroupOfFaces && nbFaces) myFacesDefaultFamilyId = REST_FACES_FAMILY;
412     if (myDoGroupOfVolumes && nbVolumes) myVolumesDefaultFamilyId = REST_VOLUMES_FAMILY;
413     if (myDoGroupOf0DElems && nb0DElements) my0DElementsDefaultFamilyId = REST_0DELEM_FAMILY;
414     if (myDoGroupOfVolumes && nbVolumes) myBallsDefaultFamilyId = REST_BALL_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          myDoGroupOf0DElems && nb0DElements,
427          myDoGroupOfBalls   && nbBalls);
428     } else {
429       aFamilies = DriverMED_Family::MakeFamilies
430         (mySubMeshes, myGroups,
431          myDoGroupOfNodes   && nbNodes,
432          myDoGroupOfEdges   && nbEdges,
433          myDoGroupOfFaces   && nbFaces,
434          myDoGroupOfVolumes && nbVolumes,
435          myDoGroupOf0DElems && nb0DElements,
436          myDoGroupOfBalls   && nbBalls);
437     }
438     //cout << " myMed->SetFamilyInfo() " << endl;
439     list<DriverMED_FamilyPtr>::iterator aFamsIter;
440     for (aFamsIter = aFamilies.begin(); aFamsIter != aFamilies.end(); aFamsIter++)
441     {
442       PFamilyInfo aFamilyInfo = (*aFamsIter)->GetFamilyInfo(myMed,aMeshInfo);
443       myMed->SetFamilyInfo(aFamilyInfo);
444     }
445
446     // Storing SMDS nodes to the MED file for the MED mesh
447     //----------------------------------------------------
448 #ifdef _EDF_NODE_IDS_
449     typedef map<TInt,TInt> TNodeIdMap;
450     TNodeIdMap aNodeIdMap;
451 #endif
452     const EModeSwitch   theMode        = eFULL_INTERLACE;
453     const ERepere       theSystem      = eCART;
454     const EBooleen      theIsElemNum   = eVRAI;
455     const EBooleen      theIsElemNames = eFAUX;
456     const EConnectivite theConnMode    = eNOD;
457
458     TInt aNbNodes = myMesh->NbNodes();
459     //cout << " myMed->CrNodeInfo() aNbNodes = " << aNbNodes << endl;
460     PNodeInfo aNodeInfo = myMed->CrNodeInfo(aMeshInfo, aNbNodes,
461                                             theMode, theSystem, theIsElemNum, theIsElemNames);
462
463     //cout << " fillElemFamilyMap( SMDSAbs_Node )" << endl;
464     // find family numbers for nodes
465     TElemFamilyMap anElemFamMap;
466     fillElemFamilyMap( anElemFamMap, aFamilies, SMDSAbs_Node );
467
468     for (TInt iNode = 0; aCoordHelperPtr->Next(); iNode++)
469     {
470       // coordinates
471       TCoordSlice aTCoordSlice = aNodeInfo->GetCoordSlice( iNode );
472       for(TInt iCoord = 0; iCoord < aSpaceDimension; iCoord++){
473         aTCoordSlice[iCoord] = aCoordHelperPtr->GetCoord(iCoord);
474       }
475       // node number
476       int aNodeID = aCoordHelperPtr->GetID();
477       aNodeInfo->SetElemNum( iNode, aNodeID );
478 #ifdef _EDF_NODE_IDS_
479       aNodeIdMap[aNodeID] = iNode+1;
480 #endif
481       // family number
482       const SMDS_MeshNode* aNode = aCoordHelperPtr->GetNode();
483       int famNum = getFamilyId( anElemFamMap, aNode, myNodesDefaultFamilyId );
484       aNodeInfo->SetFamNum( iNode, famNum );
485     }
486     anElemFamMap.Clear();
487
488     // coordinate names and units
489     for (TInt iCoord = 0; iCoord < aSpaceDimension; iCoord++) {
490       aNodeInfo->SetCoordName( iCoord, aCoordHelperPtr->GetName(iCoord));
491       aNodeInfo->SetCoordUnit( iCoord, aCoordHelperPtr->GetUnit(iCoord));
492     }
493
494     //cout << " SetNodeInfo(aNodeInfo)" << endl;
495     MESSAGE("Perform - aNodeInfo->GetNbElem() = "<<aNbNodes);
496     myMed->SetNodeInfo(aNodeInfo);
497     aNodeInfo.reset(); // free memory used for arrays
498
499
500     // Storing SMDS elements to the MED file for the MED mesh
501     //-------------------------------------------------------
502     // Write one element type at once in order to minimize memory usage (PAL19276)
503
504     const SMDS_MeshInfo& nbElemInfo = myMesh->GetMeshInfo();
505
506     // poly elements are not supported by med-2.1
507     bool polyTypesSupported = myMed->CrPolygoneInfo(aMeshInfo,eMAILLE,ePOLYGONE,0,0);
508     TInt nbPolygonNodes = 0, nbPolyhedronNodes = 0, nbPolyhedronFaces = 0;
509
510     // collect info on all geom types
511
512     list< TElemTypeData > aTElemTypeDatas;
513
514     EEntiteMaillage anEntity = eMAILLE;
515 #ifdef _ELEMENTS_BY_DIM_
516     anEntity = eNOEUD_ELEMENT;
517 #endif
518     aTElemTypeDatas.push_back(TElemTypeData(anEntity,
519                                             ePOINT1,
520                                             nbElemInfo.Nb0DElements(),
521                                             SMDSAbs_0DElement));
522 #ifdef _ELEMENTS_BY_DIM_
523     anEntity = eSTRUCT_ELEMENT;
524 #endif
525     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
526                                              eBALL,
527                                              nbElemInfo.NbBalls(),
528                                              SMDSAbs_Ball));
529 #ifdef _ELEMENTS_BY_DIM_
530     anEntity = eARETE;
531 #endif
532     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
533                                             eSEG2,
534                                             nbElemInfo.NbEdges( ORDER_LINEAR ),
535                                             SMDSAbs_Edge));
536     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
537                                             eSEG3,
538                                             nbElemInfo.NbEdges( ORDER_QUADRATIC ),
539                                             SMDSAbs_Edge));
540 #ifdef _ELEMENTS_BY_DIM_
541     anEntity = eFACE;
542 #endif
543     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
544                                             eTRIA3,
545                                             nbElemInfo.NbTriangles( ORDER_LINEAR ),
546                                             SMDSAbs_Face));
547     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
548                                             eTRIA6,
549                                             nbElemInfo.NbTriangles( ORDER_QUADRATIC ),
550                                             SMDSAbs_Face));
551     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
552                                             eQUAD4,
553                                             nbElemInfo.NbQuadrangles( ORDER_LINEAR ),
554                                             SMDSAbs_Face));
555     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
556                                              eQUAD8,
557                                              nbElemInfo.NbQuadrangles( ORDER_QUADRATIC ) -
558                                              nbElemInfo.NbBiQuadQuadrangles(),
559                                              SMDSAbs_Face));
560     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
561                                              eQUAD9,
562                                              nbElemInfo.NbBiQuadQuadrangles(),
563                                              SMDSAbs_Face));
564     if ( polyTypesSupported ) {
565       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
566                                                ePOLYGONE,
567                                                nbElemInfo.NbPolygons(),
568                                                SMDSAbs_Face));
569       // we need one more loop on poly elements to count nb of their nodes
570       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
571                                                ePOLYGONE,
572                                                nbElemInfo.NbPolygons(),
573                                                SMDSAbs_Face));
574     }
575 #ifdef _ELEMENTS_BY_DIM_
576     anEntity = eMAILLE;
577 #endif
578     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
579                                             eTETRA4,
580                                             nbElemInfo.NbTetras( ORDER_LINEAR ),
581                                             SMDSAbs_Volume));
582     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
583                                             eTETRA10,
584                                             nbElemInfo.NbTetras( ORDER_QUADRATIC ),
585                                             SMDSAbs_Volume));
586     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
587                                             ePYRA5,
588                                             nbElemInfo.NbPyramids( ORDER_LINEAR ),
589                                             SMDSAbs_Volume));
590     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
591                                             ePYRA13,
592                                             nbElemInfo.NbPyramids( ORDER_QUADRATIC ),
593                                             SMDSAbs_Volume));
594     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
595                                             ePENTA6,
596                                             nbElemInfo.NbPrisms( ORDER_LINEAR ),
597                                             SMDSAbs_Volume));
598     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
599                                              ePENTA15,
600                                              nbElemInfo.NbPrisms( ORDER_QUADRATIC ),
601                                              SMDSAbs_Volume));
602     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
603                                              eHEXA8,
604                                              nbElemInfo.NbHexas( ORDER_LINEAR ),
605                                              SMDSAbs_Volume));
606     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
607                                              eHEXA20,
608                                              nbElemInfo.NbHexas( ORDER_QUADRATIC )-
609                                              nbElemInfo.NbTriQuadHexas(),
610                                              SMDSAbs_Volume));
611     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
612                                              eHEXA27,
613                                              nbElemInfo.NbTriQuadHexas(),
614                                              SMDSAbs_Volume));
615     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
616                                              eOCTA12,
617                                              nbElemInfo.NbHexPrisms(),
618                                              SMDSAbs_Volume));
619     if ( polyTypesSupported ) {
620       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
621                                                ePOLYEDRE,
622                                                nbElemInfo.NbPolyhedrons(),
623                                                SMDSAbs_Volume));
624       // we need one more loop on poly elements to count nb of their nodes
625       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
626                                                ePOLYEDRE,
627                                                nbElemInfo.NbPolyhedrons(),
628                                                SMDSAbs_Volume));
629     }
630
631     vector< bool > isElemFamMapBuilt( SMDSAbs_NbElementTypes, false );
632
633     // loop on all geom types of elements
634
635     list< TElemTypeData >::iterator aElemTypeData = aTElemTypeDatas.begin();
636     for ( ; aElemTypeData != aTElemTypeDatas.end(); ++aElemTypeData )
637     {
638       if ( aElemTypeData->_nbElems == 0 )
639         continue;
640
641       int defaultFamilyId = 0;
642       switch ( aElemTypeData->_smdsType ) {
643       case SMDSAbs_0DElement:
644         defaultFamilyId = my0DElementsDefaultFamilyId;
645         break;
646       case SMDSAbs_Ball:
647         defaultFamilyId = myBallsDefaultFamilyId;
648         break;
649       case SMDSAbs_Edge:
650         defaultFamilyId = myEdgesDefaultFamilyId;
651         break;
652       case SMDSAbs_Face:
653         defaultFamilyId = myFacesDefaultFamilyId;
654         break;
655       case SMDSAbs_Volume:
656         defaultFamilyId = myVolumesDefaultFamilyId;
657         break;
658       default:
659         continue;
660       }
661
662       // iterator on elements of a current type
663       SMDS_ElemIteratorPtr elemIterator;
664       int iElem = 0;
665
666       // Treat POLYGONs
667       // ---------------
668       if ( aElemTypeData->_geomType == ePOLYGONE )
669       {
670         elemIterator = myMesh->elementGeomIterator( SMDSGeom_POLYGON );
671         if ( nbPolygonNodes == 0 ) {
672           // Count nb of nodes
673           while ( elemIterator->more() ) {
674             const SMDS_MeshElement* anElem = elemIterator->next();
675             nbPolygonNodes += anElem->NbNodes();
676             if ( ++iElem == aElemTypeData->_nbElems )
677               break;
678           }
679         }
680         else {
681           // Store in med file
682           PPolygoneInfo aPolygoneInfo = myMed->CrPolygoneInfo(aMeshInfo,
683                                                               aElemTypeData->_entity,
684                                                               aElemTypeData->_geomType,
685                                                               aElemTypeData->_nbElems,
686                                                               nbPolygonNodes,
687                                                               theConnMode, theIsElemNum,
688                                                               theIsElemNames);
689           TElemNum & index = *(aPolygoneInfo->myIndex.get());
690           index[0] = 1;
691
692           while ( elemIterator->more() )
693           {
694             const SMDS_MeshElement* anElem = elemIterator->next();
695             // index
696             TInt aNbNodes = anElem->NbNodes();
697             index[ iElem+1 ] = index[ iElem ] + aNbNodes;
698
699             // connectivity
700             TConnSlice aTConnSlice = aPolygoneInfo->GetConnSlice( iElem );
701             for(TInt iNode = 0; iNode < aNbNodes; iNode++) {
702               const SMDS_MeshElement* aNode = anElem->GetNode( iNode );
703 #ifdef _EDF_NODE_IDS_
704               aTConnSlice[ iNode ] = aNodeIdMap[aNode->GetID()];
705 #else
706               aTConnSlice[ iNode ] = aNode->GetID();
707 #endif
708             }
709             // element number
710             aPolygoneInfo->SetElemNum( iElem, anElem->GetID() );
711
712             // family number
713             int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
714             aPolygoneInfo->SetFamNum( iElem, famNum );
715
716             if ( ++iElem == aPolygoneInfo->GetNbElem() )
717               break;
718           }
719           myMed->SetPolygoneInfo(aPolygoneInfo);
720         }
721
722       } 
723
724       // Treat POLYEDREs
725       // ----------------
726       else if (aElemTypeData->_geomType == ePOLYEDRE )
727       {
728         elemIterator = myMesh->elementGeomIterator( SMDSGeom_POLYHEDRA );
729         
730         if ( nbPolyhedronNodes == 0 ) {
731           // Count nb of nodes
732           while ( elemIterator->more() ) {
733             const SMDS_MeshElement* anElem = elemIterator->next();
734             const SMDS_VtkVolume *aPolyedre = dynamic_cast<const SMDS_VtkVolume*>(anElem);
735             if ( !aPolyedre ) continue;
736             nbPolyhedronNodes += aPolyedre->NbNodes();
737             nbPolyhedronFaces += aPolyedre->NbFaces();
738             if ( ++iElem == aElemTypeData->_nbElems )
739               break;
740           }
741         }
742         else {
743           // Store in med file
744           PPolyedreInfo aPolyhInfo = myMed->CrPolyedreInfo(aMeshInfo,
745                                                            aElemTypeData->_entity,
746                                                            aElemTypeData->_geomType,
747                                                            aElemTypeData->_nbElems,
748                                                            nbPolyhedronFaces+1,
749                                                            nbPolyhedronNodes,
750                                                            theConnMode,
751                                                            theIsElemNum,
752                                                            theIsElemNames);
753           TElemNum & index = *(aPolyhInfo->myIndex.get());
754           TElemNum & faces = *(aPolyhInfo->myFaces.get());
755           TElemNum & conn  = *(aPolyhInfo->myConn.get());
756           index[0] = 1;
757           faces[0] = 1;
758
759           TInt iFace = 0, iNode = 0;
760           while ( elemIterator->more() )
761           {
762             const SMDS_MeshElement* anElem = elemIterator->next();
763             const SMDS_VtkVolume *aPolyedre = dynamic_cast<const SMDS_VtkVolume*>(anElem);
764             if ( !aPolyedre ) continue;
765             // index
766             TInt aNbFaces = aPolyedre->NbFaces();
767             index[ iElem+1 ] = index[ iElem ] + aNbFaces;
768
769             // face index
770             for (TInt f = 1; f <= aNbFaces; ++f, ++iFace ) {
771               int aNbFaceNodes = aPolyedre->NbFaceNodes( f );
772               faces[ iFace+1 ] = faces[ iFace ] + aNbFaceNodes;
773             }
774             // connectivity
775             SMDS_ElemIteratorPtr nodeIt = anElem->nodesIterator();
776             while ( nodeIt->more() ) {
777               const SMDS_MeshElement* aNode = nodeIt->next();
778 #ifdef _EDF_NODE_IDS_
779               conn[ iNode ] = aNodeIdMap[aNode->GetID()];
780 #else
781               conn[ iNode ] = aNode->GetID();
782 #endif
783               ++iNode;
784             }
785             // element number
786             aPolyhInfo->SetElemNum( iElem, anElem->GetID() );
787
788             // family number
789             int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
790             aPolyhInfo->SetFamNum( iElem, famNum );
791
792             if ( ++iElem == aPolyhInfo->GetNbElem() )
793               break;
794           }
795           myMed->SetPolyedreInfo(aPolyhInfo);
796         }
797       } // if (aElemTypeData->_geomType == ePOLYEDRE )
798
799       // Treat BALLs
800       // ----------------
801       else if (aElemTypeData->_geomType == eBALL )
802       {
803         // allocate data arrays
804         PBallInfo aBallInfo = myMed->CrBallInfo( aMeshInfo,
805                                                  aElemTypeData->_nbElems );
806
807         // build map of family numbers for this type
808         if ( !isElemFamMapBuilt[ aElemTypeData->_smdsType ])
809         {
810           fillElemFamilyMap( anElemFamMap, aFamilies, aElemTypeData->_smdsType );
811           isElemFamMapBuilt[ aElemTypeData->_smdsType ] = true;
812         }
813
814         elemIterator = myMesh->elementsIterator( SMDSAbs_Ball );
815         while ( elemIterator->more() )
816         {
817           const SMDS_MeshElement* anElem = elemIterator->next();
818           // connectivity
819           const SMDS_MeshElement* aNode = anElem->GetNode( 0 );
820 #ifdef _EDF_NODE_IDS_
821           (*aBallInfo->myConn)[ iElem ] = aNodeIdMap[aNode->GetID()];
822 #else
823           (*aBallInfo->myConn)[ iElem ] = aNode->GetID();
824 #endif
825           // element number
826           aBallInfo->SetElemNum( iElem, anElem->GetID() );
827
828           // diameter
829           aBallInfo->myDiameters[ iElem ] =
830             static_cast<const SMDS_BallElement*>( anElem )->GetDiameter();
831
832           // family number
833           int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
834           aBallInfo->SetFamNum( iElem, famNum );
835           ++iElem;
836         }
837         // store data in a file
838         myMed->SetBallInfo(aBallInfo);
839       }
840
841       else
842       {
843         // Treat standard types
844         // ---------------------
845
846         // allocate data arrays
847         PCellInfo aCellInfo = myMed->CrCellInfo( aMeshInfo,
848                                                  aElemTypeData->_entity,
849                                                  aElemTypeData->_geomType,
850                                                  aElemTypeData->_nbElems,
851                                                  theConnMode,
852                                                  theIsElemNum,
853                                                  theIsElemNames);
854         // build map of family numbers for this type
855         if ( !isElemFamMapBuilt[ aElemTypeData->_smdsType ])
856         {
857           //cout << " fillElemFamilyMap()" << endl;
858           fillElemFamilyMap( anElemFamMap, aFamilies, aElemTypeData->_smdsType );
859           isElemFamMapBuilt[ aElemTypeData->_smdsType ] = true;
860         }
861
862         TInt aNbNodes = MED::GetNbNodes(aElemTypeData->_geomType);
863         elemIterator = myMesh->elementsIterator( aElemTypeData->_smdsType );
864         while ( elemIterator->more() )
865         {
866           const SMDS_MeshElement* anElem = elemIterator->next();
867           if ( anElem->NbNodes() != aNbNodes || anElem->IsPoly() )
868             continue; // other geometry
869
870           // connectivity
871           TConnSlice aTConnSlice = aCellInfo->GetConnSlice( iElem );
872           for (TInt iNode = 0; iNode < aNbNodes; iNode++) {
873             const SMDS_MeshElement* aNode = anElem->GetNode( iNode );
874 #ifdef _EDF_NODE_IDS_
875             aTConnSlice[ iNode ] = aNodeIdMap[aNode->GetID()];
876 #else
877             aTConnSlice[ iNode ] = aNode->GetID();
878 #endif
879           }
880           // element number
881           aCellInfo->SetElemNum( iElem, anElem->GetID() );
882
883           // family number
884           int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
885           aCellInfo->SetFamNum( iElem, famNum );
886
887           if ( ++iElem == aCellInfo->GetNbElem() )
888             break;
889         }
890         // store data in a file
891         myMed->SetCellInfo(aCellInfo);
892       }
893
894     } // loop on geom types
895
896
897   }
898   catch(const std::exception& exc) {
899     INFOS("The following exception was caught:\n\t"<<exc.what());
900     throw;
901   }
902   catch(...) {
903     INFOS("Unknown exception was caught !!!");
904     throw;
905   }
906
907   myMeshId = -1;
908   myGroups.clear();
909   mySubMeshes.clear();
910   return aResult;
911 }