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