src/DriverMED/DriverMED_W_SMESHDS_Mesh.cxx

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