src/DriverMED/DriverMED_W_SMESHDS_Mesh.cxx

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