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