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