src/DriverMED/DriverMED_W_SMESHDS_Mesh.cxx

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