src/DriverMED/DriverMED_W_SMESHDS_Mesh.cxx

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