src/DriverMED/DriverMED_W_SMESHDS_Mesh.cxx

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