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::TIntVector::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 Driver_Mesh::Status DriverMED_W_SMESHDS_Mesh::Perform()
 347 {
 348   Status aResult = DRS_OK;
 349   try {
 350     //MESSAGE("Perform - myFile : "<<myFile);
 351
 352     // Creating the MED mesh for corresponding SMDS structure
 353     //-------------------------------------------------------
 354     string aMeshName;
 355     if (myMeshId != -1) {
 356       ostringstream aMeshNameStr;
 357       aMeshNameStr<<myMeshId;
 358       aMeshName = aMeshNameStr.str();
 359     } else {
 360       aMeshName = myMeshName;
 361     }
 362
 363     // Mesh dimension definition
 364
 365     TInt aMeshDimension = 0;
 366     if ( myMesh->NbEdges() > 0 )
 367       aMeshDimension = 1;
 368     if ( myMesh->NbFaces() > 0 )
 369       aMeshDimension = 2;
 370     if ( myMesh->NbVolumes() > 0 )
 371       aMeshDimension = 3;
 372
 373     TInt aSpaceDimension = 3;
 374     TCoordHelperPtr aCoordHelperPtr;
 375     {
 376       bool anIsXDimension = false;
 377       bool anIsYDimension = false;
 378       bool anIsZDimension = false;
 379       if ( myAutoDimension && aMeshDimension < 3 )
 380       {
 381         SMDS_NodeIteratorPtr aNodesIter = myMesh->nodesIterator();
 382         double aBounds[6] = {0.,0.,0.,0.,0.,0.};
 383         if(aNodesIter->more()){
 384           const SMDS_MeshNode* aNode = aNodesIter->next();
 385           aBounds[0] = aBounds[1] = aNode->X();
 386           aBounds[2] = aBounds[3] = aNode->Y();
 387           aBounds[4] = aBounds[5] = aNode->Z();
 388         }
 389         while(aNodesIter->more()){
 390           const SMDS_MeshNode* aNode = aNodesIter->next();
 391           aBounds[0] = min(aBounds[0],aNode->X());
 392           aBounds[1] = max(aBounds[1],aNode->X());
 393
 394           aBounds[2] = min(aBounds[2],aNode->Y());
 395           aBounds[3] = max(aBounds[3],aNode->Y());
 396
 397           aBounds[4] = min(aBounds[4],aNode->Z());
 398           aBounds[5] = max(aBounds[5],aNode->Z());
 399         }
 400
 401         double EPS = 1.0E-7;
 402         TopoDS_Shape mainShape = myMesh->ShapeToMesh();
 403         bool    hasShapeToMesh = ( myMesh->SubMeshIndices().size() > 1 );
 404         if ( !mainShape.IsNull() && hasShapeToMesh )
 405         {
 406           // define EPS by max tolerance of the mainShape (IPAL53097)
 407           TopExp_Explorer subShape;
 408           for ( subShape.Init( mainShape, TopAbs_FACE ); subShape.More(); subShape.Next() ) {
 409             EPS = Max( EPS, BRep_Tool::Tolerance( TopoDS::Face( subShape.Current() )));
 410           }
 411           for ( subShape.Init( mainShape, TopAbs_EDGE ); subShape.More(); subShape.Next() ) {
 412             EPS = Max( EPS, BRep_Tool::Tolerance( TopoDS::Edge( subShape.Current() )));
 413           }
 414           for ( subShape.Init( mainShape, TopAbs_VERTEX ); subShape.More(); subShape.Next() ) {
 415             EPS = Max( EPS, BRep_Tool::Tolerance( TopoDS::Vertex( subShape.Current() )));
 416           }
 417           EPS *= 2.;
 418         }
 419         anIsXDimension = (aBounds[1] - aBounds[0]) + abs(aBounds[1]) + abs(aBounds[0]) > EPS;
 420         anIsYDimension = (aBounds[3] - aBounds[2]) + abs(aBounds[3]) + abs(aBounds[2]) > EPS;
 421         anIsZDimension = (aBounds[5] - aBounds[4]) + abs(aBounds[5]) + abs(aBounds[4]) > EPS;
 422         if ( myZTolerance > 0 && anIsZDimension )
 423           anIsZDimension = (aBounds[5] > myZTolerance || aBounds[4] < -myZTolerance );
 424         aSpaceDimension = Max( aMeshDimension, anIsXDimension + anIsYDimension + anIsZDimension );
 425         if ( !aSpaceDimension )
 426           aSpaceDimension = 3;
 427         // PAL16857(SMESH not conform to the MED convention):
 428         if ( aSpaceDimension == 2 && anIsZDimension ) // 2D only if mesh is in XOY plane
 429           aSpaceDimension = 3;
 430         // PAL18941(a saved study with a mesh belong Z is opened and the mesh is belong X)
 431         if ( aSpaceDimension == 1 && !anIsXDimension ) {// 1D only if mesh is along OX
 432           if ( anIsYDimension ) {
 433             aSpaceDimension = 2;
 434             anIsXDimension = true;
 435           } else {
 436             aSpaceDimension = 3;
 437           }
 438         }
 439       }
 440
 441       SMDS_NodeIteratorPtr aNodesIter = myMesh->nodesIterator();
 442       switch ( aSpaceDimension ) {
 443       case 3:
 444         aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXYZGetCoord,aXYZName));
 445         break;
 446       case 2:
 447         if(anIsXDimension && anIsYDimension)
 448           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXYGetCoord,aXYName));
 449         if(anIsYDimension && anIsZDimension)
 450           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aYZGetCoord,aYZName));
 451         if(anIsXDimension && anIsZDimension)
 452           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXZGetCoord,aXZName));
 453         break;
 454       case 1:
 455         if(anIsXDimension)
 456           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aXGetCoord,aXName));
 457         if(anIsYDimension)
 458           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aYGetCoord,aYName));
 459         if(anIsZDimension)
 460           aCoordHelperPtr.reset(new TCoordHelper(aNodesIter,aZGetCoord,aZName));
 461         break;
 462       }
 463     }
 464
 465     MED::PWrapper myMed = CrWrapperW(myFile, myVersion);
 466     PMeshInfo aMeshInfo = myMed->CrMeshInfo(aMeshDimension,aSpaceDimension,aMeshName);
 467     //MESSAGE("Add - aMeshName : "<<aMeshName<<"; "<<aMeshInfo->GetName());
 468     myMed->SetMeshInfo(aMeshInfo);
 469
 470     // Storing SMDS groups and sub-meshes as med families
 471     //----------------------------------------------------
 472     int myNodesDefaultFamilyId      = 0;
 473     int my0DElementsDefaultFamilyId = 0;
 474     int myBallsDefaultFamilyId      = 0;
 475     int myEdgesDefaultFamilyId      = 0;
 476     int myFacesDefaultFamilyId      = 0;
 477     int myVolumesDefaultFamilyId    = 0;
 478     int nbNodes      = myMesh->NbNodes();
 479     int nb0DElements = myMesh->Nb0DElements();
 480     int nbBalls      = myMesh->NbBalls();
 481     int nbEdges      = myMesh->NbEdges();
 482     int nbFaces      = myMesh->NbFaces();
 483     int nbVolumes    = myMesh->NbVolumes();
 484     if (myDoGroupOfNodes)   myNodesDefaultFamilyId      = REST_NODES_FAMILY;
 485     if (myDoGroupOfEdges)   myEdgesDefaultFamilyId      = REST_EDGES_FAMILY;
 486     if (myDoGroupOfFaces)   myFacesDefaultFamilyId      = REST_FACES_FAMILY;
 487     if (myDoGroupOfVolumes) myVolumesDefaultFamilyId    = REST_VOLUMES_FAMILY;
 488     if (myDoGroupOf0DElems) my0DElementsDefaultFamilyId = REST_0DELEM_FAMILY;
 489     if (myDoGroupOfBalls)   myBallsDefaultFamilyId      = REST_BALL_FAMILY;
 490     if (myDoAllInGroups )
 491     {
 492       if (!myDoGroupOfEdges)   myEdgesDefaultFamilyId      = NIG_EDGES_FAMILY   ;
 493       if (!myDoGroupOfFaces)   myFacesDefaultFamilyId      = NIG_FACES_FAMILY   ;
 494       if (!myDoGroupOfVolumes) myVolumesDefaultFamilyId    = NIG_VOLS_FAMILY    ;
 495       if (!myDoGroupOf0DElems) my0DElementsDefaultFamilyId = NIG_0DELEM_FAMILY  ;
 496       if (!myDoGroupOfBalls)   myBallsDefaultFamilyId      = NIG_BALL_FAMILY    ;
 497     }
 498
 499     //MESSAGE("Perform - aFamilyInfo");
 500     list<DriverMED_FamilyPtr> aFamilies;
 501     if (myAllSubMeshes) {
 502       aFamilies = DriverMED_Family::MakeFamilies
 503         (myMesh->SubMeshes(), myGroups,
 504          myDoGroupOfNodes   && nbNodes,
 505          myDoGroupOfEdges   && nbEdges,
 506          myDoGroupOfFaces   && nbFaces,
 507          myDoGroupOfVolumes && nbVolumes,
 508          myDoGroupOf0DElems && nb0DElements,
 509          myDoGroupOfBalls   && nbBalls,
 510          myDoAllInGroups);
 511     }
 512     else {
 513       aFamilies = DriverMED_Family::MakeFamilies
 514         (getIterator( mySubMeshes ), myGroups,
 515          myDoGroupOfNodes   && nbNodes,
 516          myDoGroupOfEdges   && nbEdges,
 517          myDoGroupOfFaces   && nbFaces,
 518          myDoGroupOfVolumes && nbVolumes,
 519          myDoGroupOf0DElems && nb0DElements,
 520          myDoGroupOfBalls   && nbBalls,
 521          myDoAllInGroups);
 522     }
 523     list<DriverMED_FamilyPtr>::iterator aFamsIter;
 524     for (aFamsIter = aFamilies.begin(); aFamsIter != aFamilies.end(); aFamsIter++)
 525     {
 526       PFamilyInfo aFamilyInfo = (*aFamsIter)->GetFamilyInfo(myMed,aMeshInfo);
 527       myMed->SetFamilyInfo(aFamilyInfo);
 528     }
 529
 530     // Storing SMDS nodes to the MED file for the MED mesh
 531     //----------------------------------------------------
 532 #ifdef _EDF_NODE_IDS_
 533     typedef map<TInt,TInt> TNodeIdMap;
 534     TNodeIdMap aNodeIdMap;
 535 #endif
 536     const EModeSwitch   theMode        = eFULL_INTERLACE;
 537     const ERepere       theSystem      = eCART;
 538     const EBooleen      theIsElemNum   = eVRAI;
 539     const EBooleen      theIsElemNames = eFAUX;
 540     const EConnectivite theConnMode    = eNOD;
 541
 542     TInt aNbNodes = myMesh->NbNodes();
 543     PNodeInfo aNodeInfo = myMed->CrNodeInfo(aMeshInfo, aNbNodes,
 544                                             theMode, theSystem, theIsElemNum, theIsElemNames);
 545
 546     // find family numbers for nodes
 547     TElemFamilyMap anElemFamMap;
 548     fillElemFamilyMap( anElemFamMap, aFamilies, SMDSAbs_Node );
 549
 550     for (TInt iNode = 0; aCoordHelperPtr->Next(); iNode++)
 551     {
 552       // coordinates
 553       TCoordSlice aTCoordSlice = aNodeInfo->GetCoordSlice( iNode );
 554       for(TInt iCoord = 0; iCoord < aSpaceDimension; iCoord++){
 555         aTCoordSlice[iCoord] = aCoordHelperPtr->GetCoord(iCoord);
 556       }
 557       if ( aSpaceDimension == 3 &&
 558            -myZTolerance < aTCoordSlice[2] && aTCoordSlice[2] < myZTolerance )
 559         aTCoordSlice[2] = 0.;
 560
 561       // node number
 562       int aNodeID = aCoordHelperPtr->GetID();
 563       aNodeInfo->SetElemNum( iNode, aNodeID );
 564 #ifdef _EDF_NODE_IDS_
 565       aNodeIdMap.insert( aNodeIdMap.end(), make_pair( aNodeID, iNode+1 ));
 566 #endif
 567       // family number
 568       const SMDS_MeshNode* aNode = aCoordHelperPtr->GetNode();
 569       int famNum = getFamilyId( anElemFamMap, aNode, myNodesDefaultFamilyId );
 570       aNodeInfo->SetFamNum( iNode, famNum );
 571     }
 572     anElemFamMap.Clear();
 573
 574     // coordinate names and units
 575     for (TInt iCoord = 0; iCoord < aSpaceDimension; iCoord++) {
 576       aNodeInfo->SetCoordName( iCoord, aCoordHelperPtr->GetName(iCoord));
 577       aNodeInfo->SetCoordUnit( iCoord, aCoordHelperPtr->GetUnit(iCoord));
 578     }
 579
 580     //MESSAGE("Perform - aNodeInfo->GetNbElem() = "<<aNbNodes);
 581     myMed->SetNodeInfo(aNodeInfo);
 582     aNodeInfo.reset(); // free memory used for arrays
 583
 584
 585     // Storing SMDS elements to the MED file for the MED mesh
 586     //-------------------------------------------------------
 587     // Write one element type at once in order to minimize memory usage (PAL19276)
 588
 589     const SMDS_MeshInfo& nbElemInfo = myMesh->GetMeshInfo();
 590
 591     // poly elements are not supported by med-2.1
 592     bool polyTypesSupported = ( myMed->CrPolygoneInfo(aMeshInfo,eMAILLE,ePOLYGONE,0,0).get() != 0 );
 593     TInt nbPolygonNodes = 0, nbPolyhedronNodes = 0, nbPolyhedronFaces = 0;
 594
 595     // nodes on VERTEXes where 0D elements are absent
 596     std::vector<const SMDS_MeshElement*> nodesOf0D;
 597     std::vector< SMDS_ElemIteratorPtr > iterVec;
 598     SMDS_ElemIteratorPtr iterVecIter;
 599     if ( myAddODOnVertices && getNodesOfMissing0DOnVert( myMesh, nodesOf0D ))
 600     {
 601       iterVec.resize(2);
 602       iterVec[0] = myMesh->elementsIterator( SMDSAbs_0DElement );
 603       iterVec[1] = SMDS_ElemIteratorPtr
 604         ( new SMDS_ElementVectorIterator( nodesOf0D.begin(), nodesOf0D.end() ));
 605
 606       typedef SMDS_IteratorOnIterators
 607         < const SMDS_MeshElement *, std::vector< SMDS_ElemIteratorPtr > > TItIterator;
 608       iterVecIter = SMDS_ElemIteratorPtr( new TItIterator( iterVec ));
 609     }
 610
 611     // collect info on all geom types
 612
 613     list< TElemTypeData > aTElemTypeDatas;
 614
 615     EEntiteMaillage anEntity = eMAILLE;
 616 #ifdef _ELEMENTS_BY_DIM_
 617     anEntity = eNOEUD_ELEMENT;
 618 #endif
 619     aTElemTypeDatas.push_back(TElemTypeData(anEntity,
 620                                             ePOINT1,
 621                                             nbElemInfo.Nb0DElements() + nodesOf0D.size(),
 622                                             SMDSAbs_0DElement));
 623 #ifdef _ELEMENTS_BY_DIM_
 624     anEntity = eSTRUCT_ELEMENT;
 625 #endif
 626     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 627                                              eBALL,
 628                                              nbElemInfo.NbBalls(),
 629                                              SMDSAbs_Ball));
 630 #ifdef _ELEMENTS_BY_DIM_
 631     anEntity = eARETE;
 632 #endif
 633     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 634                                              eSEG2,
 635                                              nbElemInfo.NbEdges( ORDER_LINEAR ),
 636                                              SMDSAbs_Edge));
 637     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 638                                              eSEG3,
 639                                              nbElemInfo.NbEdges( ORDER_QUADRATIC ),
 640                                              SMDSAbs_Edge));
 641 #ifdef _ELEMENTS_BY_DIM_
 642     anEntity = eFACE;
 643 #endif
 644     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 645                                              eTRIA3,
 646                                              nbElemInfo.NbTriangles( ORDER_LINEAR ),
 647                                              SMDSAbs_Face));
 648     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 649                                              eTRIA6,
 650                                              nbElemInfo.NbTriangles( ORDER_QUADRATIC ) -
 651                                              nbElemInfo.NbBiQuadTriangles(),
 652                                              SMDSAbs_Face));
 653     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 654                                              eTRIA7,
 655                                              nbElemInfo.NbBiQuadTriangles(),
 656                                              SMDSAbs_Face));
 657     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 658                                              eQUAD4,
 659                                              nbElemInfo.NbQuadrangles( ORDER_LINEAR ),
 660                                              SMDSAbs_Face));
 661     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 662                                              eQUAD8,
 663                                              nbElemInfo.NbQuadrangles( ORDER_QUADRATIC ) -
 664                                              nbElemInfo.NbBiQuadQuadrangles(),
 665                                              SMDSAbs_Face));
 666     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 667                                              eQUAD9,
 668                                              nbElemInfo.NbBiQuadQuadrangles(),
 669                                              SMDSAbs_Face));
 670     if ( polyTypesSupported ) {
 671       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 672                                                ePOLYGONE,
 673                                                nbElemInfo.NbPolygons( ORDER_LINEAR ),
 674                                                SMDSAbs_Face));
 675       // we need one more loop on poly elements to count nb of their nodes
 676       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 677                                                ePOLYGONE,
 678                                                nbElemInfo.NbPolygons( ORDER_LINEAR ),
 679                                                SMDSAbs_Face));
 680       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 681                                                ePOLYGON2,
 682                                                nbElemInfo.NbPolygons( ORDER_QUADRATIC ),
 683                                                SMDSAbs_Face));
 684       // we need one more loop on QUAD poly elements to count nb of their nodes
 685       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 686                                                ePOLYGON2,
 687                                                nbElemInfo.NbPolygons( ORDER_QUADRATIC ),
 688                                                SMDSAbs_Face));
 689     }
 690 #ifdef _ELEMENTS_BY_DIM_
 691     anEntity = eMAILLE;
 692 #endif
 693     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 694                                              eTETRA4,
 695                                              nbElemInfo.NbTetras( ORDER_LINEAR ),
 696                                              SMDSAbs_Volume));
 697     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 698                                              eTETRA10,
 699                                              nbElemInfo.NbTetras( ORDER_QUADRATIC ),
 700                                              SMDSAbs_Volume));
 701     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 702                                              ePYRA5,
 703                                              nbElemInfo.NbPyramids( ORDER_LINEAR ),
 704                                              SMDSAbs_Volume));
 705     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 706                                              ePYRA13,
 707                                              nbElemInfo.NbPyramids( ORDER_QUADRATIC ),
 708                                              SMDSAbs_Volume));
 709     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 710                                              ePENTA6,
 711                                              nbElemInfo.NbPrisms( ORDER_LINEAR ),
 712                                              SMDSAbs_Volume));
 713     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 714                                              ePENTA15,
 715                                              nbElemInfo.NbQuadPrisms(),
 716                                              SMDSAbs_Volume));
 717     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 718                                              ePENTA18,
 719                                              nbElemInfo.NbBiQuadPrisms(),
 720                                              SMDSAbs_Volume));
 721     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 722                                              eHEXA8,
 723                                              nbElemInfo.NbHexas( ORDER_LINEAR ),
 724                                              SMDSAbs_Volume));
 725     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 726                                              eHEXA20,
 727                                              nbElemInfo.NbHexas( ORDER_QUADRATIC )-
 728                                              nbElemInfo.NbTriQuadHexas(),
 729                                              SMDSAbs_Volume));
 730     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 731                                              eHEXA27,
 732                                              nbElemInfo.NbTriQuadHexas(),
 733                                              SMDSAbs_Volume));
 734     aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 735                                              eOCTA12,
 736                                              nbElemInfo.NbHexPrisms(),
 737                                              SMDSAbs_Volume));
 738     if ( polyTypesSupported ) {
 739       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 740                                                ePOLYEDRE,
 741                                                nbElemInfo.NbPolyhedrons(),
 742                                                SMDSAbs_Volume));
 743       // we need one more loop on poly elements to count nb of their nodes
 744       aTElemTypeDatas.push_back( TElemTypeData(anEntity,
 745                                                ePOLYEDRE,
 746                                                nbElemInfo.NbPolyhedrons(),
 747                                                SMDSAbs_Volume));
 748     }
 749
 750     vector< bool > isElemFamMapBuilt( SMDSAbs_NbElementTypes, false );
 751
 752     // loop on all geom types of elements
 753
 754     list< TElemTypeData >::iterator aElemTypeData = aTElemTypeDatas.begin();
 755     for ( ; aElemTypeData != aTElemTypeDatas.end(); ++aElemTypeData )
 756     {
 757       if ( aElemTypeData->_nbElems == 0 )
 758         continue;
 759
 760       int defaultFamilyId = 0;
 761       switch ( aElemTypeData->_smdsType ) {
 762       case SMDSAbs_0DElement: defaultFamilyId = my0DElementsDefaultFamilyId; break;
 763       case SMDSAbs_Ball:      defaultFamilyId = myBallsDefaultFamilyId;      break;
 764       case SMDSAbs_Edge:      defaultFamilyId = myEdgesDefaultFamilyId;      break;
 765       case SMDSAbs_Face:      defaultFamilyId = myFacesDefaultFamilyId;      break;
 766       case SMDSAbs_Volume:    defaultFamilyId = myVolumesDefaultFamilyId;    break;
 767       default:
 768         continue;
 769       }
 770
 771       // build map of family numbers for this type
 772       if ( !isElemFamMapBuilt[ aElemTypeData->_smdsType ])
 773       {
 774         fillElemFamilyMap( anElemFamMap, aFamilies, aElemTypeData->_smdsType );
 775         isElemFamMapBuilt[ aElemTypeData->_smdsType ] = true;
 776       }
 777
 778       // iterator on elements of a current type
 779       SMDS_ElemIteratorPtr elemIterator;
 780       int iElem = 0;
 781
 782       // Treat POLYGONs
 783       // ---------------
 784       if ( aElemTypeData->_geomType == ePOLYGONE ||
 785            aElemTypeData->_geomType == ePOLYGON2 )
 786       {
 787         if ( aElemTypeData->_geomType == ePOLYGONE )
 788           elemIterator = myMesh->elementEntityIterator( SMDSEntity_Polygon );
 789         else
 790           elemIterator = myMesh->elementEntityIterator( SMDSEntity_Quad_Polygon );
 791
 792         if ( nbPolygonNodes == 0 ) {
 793           // Count nb of nodes
 794           while ( elemIterator->more() ) {
 795             const SMDS_MeshElement* anElem = elemIterator->next();
 796             nbPolygonNodes += anElem->NbNodes();
 797             if ( ++iElem == aElemTypeData->_nbElems )
 798               break;
 799           }
 800         }
 801         else {
 802           // Store in med file
 803           PPolygoneInfo aPolygoneInfo = myMed->CrPolygoneInfo(aMeshInfo,
 804                                                               aElemTypeData->_entity,
 805                                                               aElemTypeData->_geomType,
 806                                                               aElemTypeData->_nbElems,
 807                                                               nbPolygonNodes,
 808                                                               theConnMode, theIsElemNum,
 809                                                               theIsElemNames);
 810           TElemNum & index = *(aPolygoneInfo->myIndex.get());
 811           index[0] = 1;
 812
 813           while ( elemIterator->more() )
 814           {
 815             const SMDS_MeshElement* anElem = elemIterator->next();
 816             // index
 817             TInt aNbNodes = anElem->NbNodes();
 818             index[ iElem+1 ] = index[ iElem ] + aNbNodes;
 819
 820             // connectivity
 821             TConnSlice aTConnSlice = aPolygoneInfo->GetConnSlice( iElem );
 822             for(TInt iNode = 0; iNode < aNbNodes; iNode++) {
 823               const SMDS_MeshElement* aNode = anElem->GetNode( iNode );
 824 #ifdef _EDF_NODE_IDS_
 825               aTConnSlice[ iNode ] = aNodeIdMap[aNode->GetID()];
 826 #else
 827               aTConnSlice[ iNode ] = aNode->GetID();
 828 #endif
 829             }
 830             // element number
 831             aPolygoneInfo->SetElemNum( iElem, anElem->GetID() );
 832
 833             // family number
 834             int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
 835             aPolygoneInfo->SetFamNum( iElem, famNum );
 836
 837             if ( ++iElem == aPolygoneInfo->GetNbElem() )
 838               break;
 839           }
 840           myMed->SetPolygoneInfo(aPolygoneInfo);
 841
 842           nbPolygonNodes = 0; // to treat next polygon type
 843         }
 844       }
 845
 846       // Treat POLYEDREs
 847       // ----------------
 848       else if ( aElemTypeData->_geomType == ePOLYEDRE )
 849       {
 850         elemIterator = myMesh->elementGeomIterator( SMDSGeom_POLYHEDRA );
 851
 852         if ( nbPolyhedronNodes == 0 ) {
 853           // Count nb of nodes
 854           while ( elemIterator->more() ) {
 855             const SMDS_MeshElement*  anElem = elemIterator->next();
 856             nbPolyhedronNodes += anElem->NbNodes();
 857             nbPolyhedronFaces += anElem->NbFaces();
 858             if ( ++iElem == aElemTypeData->_nbElems )
 859               break;
 860           }
 861         }
 862         else {
 863           // Store in med file
 864           PPolyedreInfo aPolyhInfo = myMed->CrPolyedreInfo(aMeshInfo,
 865                                                            aElemTypeData->_entity,
 866                                                            aElemTypeData->_geomType,
 867                                                            aElemTypeData->_nbElems,
 868                                                            nbPolyhedronFaces+1,
 869                                                            nbPolyhedronNodes,
 870                                                            theConnMode,
 871                                                            theIsElemNum,
 872                                                            theIsElemNames);
 873           TElemNum & index = *(aPolyhInfo->myIndex.get());
 874           TElemNum & faces = *(aPolyhInfo->myFaces.get());
 875           TElemNum & conn  = *(aPolyhInfo->myConn.get());
 876           index[0] = 1;
 877           faces[0] = 1;
 878
 879           TInt iFace = 0, iNode = 0;
 880           while ( elemIterator->more() )
 881           {
 882             const SMDS_MeshElement* anElem = elemIterator->next();
 883             const SMDS_MeshVolume *aPolyedre = myMesh->DownCast< SMDS_MeshVolume >( anElem );
 884             if ( !aPolyedre ) continue;
 885             // index
 886             TInt aNbFaces = aPolyedre->NbFaces();
 887             index[ iElem+1 ] = index[ iElem ] + aNbFaces;
 888
 889             // face index
 890             for (TInt f = 1; f <= aNbFaces; ++f, ++iFace ) {
 891               int aNbFaceNodes = aPolyedre->NbFaceNodes( f );
 892               faces[ iFace+1 ] = faces[ iFace ] + aNbFaceNodes;
 893             }
 894             // connectivity
 895             SMDS_ElemIteratorPtr nodeIt = anElem->nodesIterator();
 896             while ( nodeIt->more() ) {
 897               const SMDS_MeshElement* aNode = nodeIt->next();
 898 #ifdef _EDF_NODE_IDS_
 899               conn[ iNode ] = aNodeIdMap[aNode->GetID()];
 900 #else
 901               conn[ iNode ] = aNode->GetID();
 902 #endif
 903               ++iNode;
 904             }
 905             // element number
 906             aPolyhInfo->SetElemNum( iElem, anElem->GetID() );
 907
 908             // family number
 909             int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
 910             aPolyhInfo->SetFamNum( iElem, famNum );
 911
 912             if ( ++iElem == aPolyhInfo->GetNbElem() )
 913               break;
 914           }
 915           myMed->SetPolyedreInfo(aPolyhInfo);
 916         }
 917       } // if (aElemTypeData->_geomType == ePOLYEDRE )
 918
 919       // Treat BALLs
 920       // ----------------
 921       else if (aElemTypeData->_geomType == eBALL )
 922       {
 923         // allocate data arrays
 924         PBallInfo aBallInfo = myMed->CrBallInfo( aMeshInfo, aElemTypeData->_nbElems );
 925
 926         elemIterator = myMesh->elementsIterator( SMDSAbs_Ball );
 927         while ( elemIterator->more() )
 928         {
 929           const SMDS_MeshElement* anElem = elemIterator->next();
 930           // connectivity
 931           const SMDS_MeshElement* aNode = anElem->GetNode( 0 );
 932 #ifdef _EDF_NODE_IDS_
 933           (*aBallInfo->myConn)[ iElem ] = aNodeIdMap[aNode->GetID()];
 934 #else
 935           (*aBallInfo->myConn)[ iElem ] = aNode->GetID();
 936 #endif
 937           // element number
 938           aBallInfo->SetElemNum( iElem, anElem->GetID() );
 939
 940           // diameter
 941           aBallInfo->myDiameters[ iElem ] =
 942             static_cast<const SMDS_BallElement*>( anElem )->GetDiameter();
 943
 944           // family number
 945           int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
 946           aBallInfo->SetFamNum( iElem, famNum );
 947           ++iElem;
 948         }
 949         // store data in a file
 950         myMed->SetBallInfo(aBallInfo);
 951       }
 952
 953       else
 954       {
 955         // Treat standard types
 956         // ---------------------
 957         // allocate data arrays
 958         PCellInfo aCellInfo = myMed->CrCellInfo( aMeshInfo,
 959                                                  aElemTypeData->_entity,
 960                                                  aElemTypeData->_geomType,
 961                                                  aElemTypeData->_nbElems,
 962                                                  theConnMode,
 963                                                  theIsElemNum,
 964                                                  theIsElemNames);
 965
 966         TInt aNbNodes = MED::GetNbNodes(aElemTypeData->_geomType);
 967         elemIterator = myMesh->elementsIterator( aElemTypeData->_smdsType );
 968         if ( aElemTypeData->_smdsType == SMDSAbs_0DElement && ! nodesOf0D.empty() )
 969           elemIterator = iterVecIter;
 970         while ( elemIterator->more() )
 971         {
 972           const SMDS_MeshElement* anElem = elemIterator->next();
 973           if ( anElem->NbNodes() != aNbNodes || anElem->IsPoly() )
 974             continue; // other geometry
 975
 976           // connectivity
 977           TConnSlice aTConnSlice = aCellInfo->GetConnSlice( iElem );
 978           for (TInt iNode = 0; iNode < aNbNodes; iNode++) {
 979             const SMDS_MeshElement* aNode = anElem->GetNode( iNode );
 980 #ifdef _EDF_NODE_IDS_
 981             aTConnSlice[ iNode ] = aNodeIdMap[aNode->GetID()];
 982 #else
 983             aTConnSlice[ iNode ] = aNode->GetID();
 984 #endif
 985           }
 986           // element number
 987           aCellInfo->SetElemNum( iElem, anElem->GetID() );
 988
 989           // family number
 990           int famNum = getFamilyId( anElemFamMap, anElem, defaultFamilyId );
 991           aCellInfo->SetFamNum( iElem, famNum );
 992
 993           if ( ++iElem == aCellInfo->GetNbElem() )
 994             break;
 995         }
 996         // store data in a file
 997         myMed->SetCellInfo(aCellInfo);
 998       }
 999     } // loop on geom types
1000   }
1001   catch(const std::exception& exc) {
1002     INFOS("The following exception was caught:\n\t"<<exc.what());
1003     throw;
1004   }
1005   catch(...) {
1006     INFOS("Unknown exception was caught !!!");
1007     throw;
1008   }
1009
1010   myMeshId = -1;
1011   myGroups.clear();
1012   mySubMeshes.clear();
1013   return aResult;
1014 }
1015
1016 //================================================================================
1017 /*!
1018  * \brief Returns nodes on VERTEXes where 0D elements are absent
1019  */
1020 //================================================================================
1021
1022 bool DriverMED_W_SMESHDS_Mesh::
1023 getNodesOfMissing0DOnVert(SMESHDS_Mesh*                         meshDS,
1024                           std::vector<const SMDS_MeshElement*>& nodes)
1025 {
1026   nodes.clear();
1027   for ( int i = 1; i <= meshDS->MaxShapeIndex(); ++i )
1028   {
1029     if ( meshDS->IndexToShape( i ).ShapeType() != TopAbs_VERTEX )
1030       continue;
1031     if ( SMESHDS_SubMesh* sm = meshDS->MeshElements(i) ) {
1032       SMDS_NodeIteratorPtr nIt= sm->GetNodes();
1033       while (nIt->more())
1034       {
1035         const SMDS_MeshNode* n = nIt->next();
1036         if ( n->NbInverseElements( SMDSAbs_0DElement ) == 0 )
1037           nodes.push_back( n );
1038       }
1039     }
1040   }
1041   return !nodes.empty();
1042 }