src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

   1 // Copyright (C) 2004-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 //=============================================================================
  21 // File      : GHS3DPlugin_GHS3D.cxx
  22 // Created   :
  23 // Author    : Edward AGAPOV, modified by Lioka RAZAFINDRAZAKA (CEA) 09/02/2007
  24 // Project   : SALOME
  25 //=============================================================================
  26 //
  27 #include "GHS3DPlugin_GHS3D.hxx"
  28 #include "GHS3DPlugin_Hypothesis.hxx"
  29
  30 #include <Basics_Utils.hxx>
  31
  32 //#include <SMESH_Gen.hxx>
  33 #include <SMESH_Client.hxx>
  34 #include <SMESH_Mesh.hxx>
  35 #include <SMESH_Comment.hxx>
  36 #include <SMESH_MesherHelper.hxx>
  37 #include <SMESH_MeshEditor.hxx>
  38 #include <SMESH_OctreeNode.hxx>
  39 #include <SMESH_Group.hxx>
  40 #include <SMESH_subMeshEventListener.hxx>
  41 #include <SMESH_HypoFilter.hxx>
  42 #include <SMESH_MeshAlgos.hxx>
  43
  44 #include <SMDS_MeshElement.hxx>
  45 #include <SMDS_MeshNode.hxx>
  46 #include <SMDS_FaceOfNodes.hxx>
  47 #include <SMDS_VolumeOfNodes.hxx>
  48
  49 #include <SMESHDS_Group.hxx>
  50
  51 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  52 #include <StdMeshers_ViscousLayers.hxx>
  53
  54 #include <BRepAdaptor_Surface.hxx>
  55 #include <BRepBndLib.hxx>
  56 #include <BRepBuilderAPI_MakeVertex.hxx>
  57 #include <BRepClass3d_SolidClassifier.hxx>
  58 #include <BRepExtrema_DistShapeShape.hxx>
  59 #include <BRepGProp.hxx>
  60 #include <BRepTools.hxx>
  61 #include <BRep_Tool.hxx>
  62 #include <Bnd_Box.hxx>
  63 #include <GProp_GProps.hxx>
  64 #include <GeomAPI_ProjectPointOnSurf.hxx>
  65 #include <OSD_File.hxx>
  66 #include <Precision.hxx>
  67 #include <Quantity_Parameter.hxx>
  68 #include <Standard_ErrorHandler.hxx>
  69 #include <Standard_Failure.hxx>
  70 #include <Standard_ProgramError.hxx>
  71 #include <TopExp.hxx>
  72 #include <TopExp_Explorer.hxx>
  73 #include <TopTools_IndexedMapOfShape.hxx>
  74 #include <TopTools_ListIteratorOfListOfShape.hxx>
  75 #include <TopTools_MapOfShape.hxx>
  76 #include <TopoDS.hxx>
  77 #include <TopoDS_Shape.hxx>
  78 #include <TopoDS_Solid.hxx>
  79
  80 #include <utilities.h>
  81
  82 #ifdef WIN32
  83 #include <io.h>
  84 #else
  85 #include <sys/sysinfo.h>
  86 #endif
  87 #include <algorithm>
  88
  89 //#include <Standard_Stream.hxx>
  90
  91
  92 #define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
  93
  94 #ifdef _DEBUG_
  95 #define DUMP(txt) \
  96 //  std::cout << txt
  97 #else
  98 #define DUMP(txt)
  99 #endif
 100
 101 extern "C"
 102 {
 103 #ifndef WNT
 104 #include <unistd.h>
 105 #include <sys/mman.h>
 106 #endif
 107 #include <sys/stat.h>
 108 #include <fcntl.h>
 109 }
 110
 111 #define HOLE_ID -1
 112
 113 typedef const list<const SMDS_MeshFace*> TTriaList;
 114
 115 static const char theDomainGroupNamePrefix[] = "Domain_";
 116
 117 static void removeFile( const TCollection_AsciiString& fileName )
 118 {
 119   try {
 120     OSD_File( fileName ).Remove();
 121   }
 122   catch ( Standard_ProgramError ) {
 123     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
 124   }
 125 }
 126
 127 //=============================================================================
 128 /*!
 129  *
 130  */
 131 //=============================================================================
 132
 133 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
 134   : SMESH_3D_Algo(hypId, studyId, gen)
 135 {
 136   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
 137   _name = Name();
 138   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
 139   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
 140   _iShape=0;
 141   _nbShape=0;
 142   _compatibleHypothesis.push_back( GHS3DPlugin_Hypothesis::GetHypType());
 143   _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() );
 144   _requireShape = false; // can work without shape_studyId
 145
 146   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 147   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 148   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 149
 150   MESSAGE("studyid = " << _studyId);
 151
 152   myStudy = NULL;
 153   myStudy = aStudyMgr->GetStudyByID(_studyId);
 154   if (myStudy)
 155     MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
 156
 157   _compute_canceled = false;
 158 }
 159
 160 //=============================================================================
 161 /*!
 162  *
 163  */
 164 //=============================================================================
 165
 166 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
 167 {
 168   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
 169 }
 170
 171 //=============================================================================
 172 /*!
 173  *
 174  */
 175 //=============================================================================
 176
 177 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&         aMesh,
 178                                           const TopoDS_Shape& aShape,
 179                                           Hypothesis_Status&  aStatus )
 180 {
 181   aStatus = SMESH_Hypothesis::HYP_OK;
 182
 183   _hyp = 0;
 184   _viscousLayersHyp = 0;
 185   _keepFiles = false;
 186   _removeLogOnSuccess = true;
 187   _logInStandardOutput = false;
 188
 189   const list <const SMESHDS_Hypothesis * >& hyps =
 190     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 191   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 192   for ( ; h != hyps.end(); ++h )
 193   {
 194     if ( !_hyp )
 195       _hyp = dynamic_cast< const GHS3DPlugin_Hypothesis*> ( *h );
 196     if ( !_viscousLayersHyp )
 197       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 198   }
 199   if ( _hyp )
 200   {
 201     _keepFiles = _hyp->GetKeepFiles();
 202     _removeLogOnSuccess = _hyp->GetRemoveLogOnSuccess();
 203     _logInStandardOutput = _hyp->GetStandardOutputLog();
 204   }
 205
 206   return true;
 207 }
 208
 209
 210 //=======================================================================
 211 //function : entryToShape
 212 //purpose  :
 213 //=======================================================================
 214
 215 TopoDS_Shape GHS3DPlugin_GHS3D::entryToShape(std::string entry)
 216 {
 217   MESSAGE("GHS3DPlugin_GHS3D::entryToShape "<<entry );
 218   GEOM::GEOM_Object_var aGeomObj;
 219   TopoDS_Shape S = TopoDS_Shape();
 220   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 221   if (!aSObj->_is_nil() ) {
 222     CORBA::Object_var obj = aSObj->GetObject();
 223     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 224     aSObj->UnRegister();
 225   }
 226   if ( !aGeomObj->_is_nil() )
 227     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 228   return S;
 229 }
 230
 231 //=======================================================================
 232 //function : findShape
 233 //purpose  :
 234 //=======================================================================
 235
 236 static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[],
 237                               TopoDS_Shape        aShape,
 238                               const TopoDS_Shape  shape[],
 239                               double**            box,
 240                               const int           nShape,
 241                               TopAbs_State *      state = 0)
 242 {
 243   gp_XYZ aPnt(0,0,0);
 244   int j, iShape, nbNode = 4;
 245
 246   for ( j=0; j<nbNode; j++ ) {
 247     gp_XYZ p ( aNode[j]->X(), aNode[j]->Y(), aNode[j]->Z() );
 248     if ( aNode[j]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) {
 249       aPnt = p;
 250       break;
 251     }
 252     aPnt += p / nbNode;
 253   }
 254
 255   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 256   if (state) *state = SC.State();
 257   if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) {
 258     for (iShape = 0; iShape < nShape; iShape++) {
 259       aShape = shape[iShape];
 260       if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() ||
 261               aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() ||
 262               aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) {
 263         BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 264         if (state) *state = SC.State();
 265         if (SC.State() == TopAbs_IN)
 266           break;
 267       }
 268     }
 269   }
 270   return aShape;
 271 }
 272
 273 //=======================================================================
 274 //function : readMapIntLine
 275 //purpose  :
 276 //=======================================================================
 277
 278 static char* readMapIntLine(char* ptr, int tab[]) {
 279   long int intVal;
 280   std::cout << std::endl;
 281
 282   for ( int i=0; i<17; i++ ) {
 283     intVal = strtol(ptr, &ptr, 10);
 284     if ( i < 3 )
 285       tab[i] = intVal;
 286   }
 287   return ptr;
 288 }
 289
 290 //================================================================================
 291 /*!
 292  * \brief returns true if a triangle defined by the nodes is a temporary face on a
 293  * side facet of pyramid and defines sub-domian inside the pyramid
 294  */
 295 //================================================================================
 296
 297 static bool isTmpFace(const SMDS_MeshNode* node1,
 298                       const SMDS_MeshNode* node2,
 299                       const SMDS_MeshNode* node3)
 300 {
 301   // find a pyramid sharing the 3 nodes
 302   //const SMDS_MeshElement* pyram = 0;
 303   SMDS_ElemIteratorPtr vIt1 = node1->GetInverseElementIterator(SMDSAbs_Volume);
 304   while ( vIt1->more() )
 305   {
 306     const SMDS_MeshElement* pyram = vIt1->next();
 307     if ( pyram->NbCornerNodes() != 5 ) continue;
 308     int i2, i3;
 309     if ( (i2 = pyram->GetNodeIndex( node2 )) >= 0 &&
 310          (i3 = pyram->GetNodeIndex( node3 )) >= 0 )
 311     {
 312       // Triangle defines sub-domian inside the pyramid if it's
 313       // normal points out of the pyram
 314
 315       // make i2 and i3 hold indices of base nodes of the pyram while
 316       // keeping the nodes order in the triangle
 317       const int iApex = 4;
 318       if ( i2 == iApex )
 319         i2 = i3, i3 = pyram->GetNodeIndex( node1 );
 320       else if ( i3 == iApex )
 321         i3 = i2, i2 = pyram->GetNodeIndex( node1 );
 322
 323       int i3base = (i2+1) % 4; // next index after i2 within the pyramid base
 324       return ( i3base != i3 );
 325     }
 326   }
 327   return false;
 328 }
 329
 330 //=======================================================================
 331 //function : findShapeID
 332 //purpose  : find the solid corresponding to GHS3D sub-domain following
 333 //           the technique proposed in GHS3D manual (available within
 334 //           ghs3d installation) in chapter "B.4 Subdomain (sub-region) assignment".
 335 //           In brief: normal of the triangle defined by the given nodes
 336 //           points out of the domain it is associated to
 337 //=======================================================================
 338
 339 static int findShapeID(SMESH_Mesh&          mesh,
 340                        const SMDS_MeshNode* node1,
 341                        const SMDS_MeshNode* node2,
 342                        const SMDS_MeshNode* node3,
 343                        const bool           toMeshHoles)
 344 {
 345   const int invalidID = 0;
 346   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 347
 348   // face the nodes belong to
 349   const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
 350   if ( !face )
 351     return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
 352 #ifdef _DEBUG_
 353   std::cout << "bnd face " << face->GetID() << " - ";
 354 #endif
 355   // geom face the face assigned to
 356   SMESH_MeshEditor editor(&mesh);
 357   int geomFaceID = editor.FindShape( face );
 358   if ( !geomFaceID )
 359     return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
 360   TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
 361   if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
 362     return invalidID;
 363   TopoDS_Face geomFace = TopoDS::Face( shape );
 364
 365   // solids bounded by geom face
 366   TopTools_IndexedMapOfShape solids, shells;
 367   TopTools_ListIteratorOfListOfShape ansIt = mesh.GetAncestors(geomFace);
 368   for ( ; ansIt.More(); ansIt.Next() ) {
 369     switch ( ansIt.Value().ShapeType() ) {
 370     case TopAbs_SOLID:
 371       solids.Add( ansIt.Value() ); break;
 372     case TopAbs_SHELL:
 373       shells.Add( ansIt.Value() ); break;
 374     default:;
 375     }
 376   }
 377   // analyse found solids
 378   if ( solids.Extent() == 0 || shells.Extent() == 0)
 379     return invalidID;
 380
 381   const TopoDS_Solid& solid1 = TopoDS::Solid( solids(1) );
 382   if ( solids.Extent() == 1 )
 383   {
 384     if ( toMeshHoles )
 385       return meshDS->ShapeToIndex( solid1 );
 386
 387     //////////// UNCOMMENT AS SOON AS
 388     //////////// http://tracker.dev.opencascade.org/view.php?id=23129
 389     //////////// IS SOLVED
 390     // - Are we at a hole boundary face?
 391     // if ( shells(1).IsSame( BRepTools::OuterShell( solid1 )) )
 392     // { // - No, but maybe a hole is bound by two shapes? Does shells(1) touches another shell?
 393     //   bool touch = false;
 394     //   TopExp_Explorer eExp( shells(1), TopAbs_EDGE );
 395     //   // check if any edge of shells(1) belongs to another shell
 396     //   for ( ; eExp.More() && !touch; eExp.Next() ) {
 397     //     ansIt = mesh.GetAncestors( eExp.Current() );
 398     //     for ( ; ansIt.More() && !touch; ansIt.Next() ) {
 399     //       if ( ansIt.Value().ShapeType() == TopAbs_SHELL )
 400     //         touch = ( !ansIt.Value().IsSame( shells(1) ));
 401     //     }
 402     //   }
 403     //   if (!touch)
 404     //     return meshDS->ShapeToIndex( solid1 );
 405     // }
 406   }
 407   // find orientation of geom face within the first solid
 408   TopExp_Explorer fExp( solid1, TopAbs_FACE );
 409   for ( ; fExp.More(); fExp.Next() )
 410     if ( geomFace.IsSame( fExp.Current() )) {
 411       geomFace = TopoDS::Face( fExp.Current() );
 412       break;
 413     }
 414   if ( !fExp.More() )
 415     return invalidID; // face not found
 416
 417   // normale to triangle
 418   gp_Pnt node1Pnt ( node1->X(), node1->Y(), node1->Z() );
 419   gp_Pnt node2Pnt ( node2->X(), node2->Y(), node2->Z() );
 420   gp_Pnt node3Pnt ( node3->X(), node3->Y(), node3->Z() );
 421   gp_Vec vec12( node1Pnt, node2Pnt );
 422   gp_Vec vec13( node1Pnt, node3Pnt );
 423   gp_Vec meshNormal = vec12 ^ vec13;
 424   if ( meshNormal.SquareMagnitude() < DBL_MIN )
 425     return invalidID;
 426
 427   // get normale to geomFace at any node
 428   bool geomNormalOK = false;
 429   gp_Vec geomNormal;
 430   const SMDS_MeshNode* nodes[3] = { node1, node2, node3 };
 431   SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace );
 432   for ( int i = 0; !geomNormalOK && i < 3; ++i )
 433   {
 434     // find UV of i-th node on geomFace
 435     const SMDS_MeshNode* nNotOnSeamEdge = 0;
 436     if ( helper.IsSeamShape( nodes[i]->getshapeId() )) {
 437       if ( helper.IsSeamShape( nodes[(i+1)%3]->getshapeId() ))
 438         nNotOnSeamEdge = nodes[(i+2)%3];
 439       else
 440         nNotOnSeamEdge = nodes[(i+1)%3];
 441     }
 442     bool uvOK;
 443     gp_XY uv = helper.GetNodeUV( geomFace, nodes[i], nNotOnSeamEdge, &uvOK );
 444     // check that uv is correct
 445     if (uvOK) {
 446       double tol = 1e-6;
 447       TopoDS_Shape nodeShape = helper.GetSubShapeByNode( nodes[i], meshDS );
 448       if ( !nodeShape.IsNull() )
 449         switch ( nodeShape.ShapeType() )
 450         {
 451         case TopAbs_FACE:   tol = BRep_Tool::Tolerance( TopoDS::Face( nodeShape )); break;
 452         case TopAbs_EDGE:   tol = BRep_Tool::Tolerance( TopoDS::Edge( nodeShape )); break;
 453         case TopAbs_VERTEX: tol = BRep_Tool::Tolerance( TopoDS::Vertex( nodeShape )); break;
 454         default:;
 455         }
 456       gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
 457       BRepAdaptor_Surface surface( geomFace );
 458       uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
 459       if ( uvOK ) {
 460         // normale to geomFace at UV
 461         gp_Vec du, dv;
 462         surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
 463         geomNormal = du ^ dv;
 464         if ( geomFace.Orientation() == TopAbs_REVERSED )
 465           geomNormal.Reverse();
 466         geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
 467       }
 468     }
 469   }
 470   if ( !geomNormalOK)
 471     return invalidID;
 472
 473   // compare normals
 474   bool isReverse = ( meshNormal * geomNormal ) < 0;
 475   if ( !isReverse )
 476     return meshDS->ShapeToIndex( solid1 );
 477
 478   if ( solids.Extent() == 1 )
 479     return HOLE_ID; // we are inside a hole
 480   else
 481     return meshDS->ShapeToIndex( solids(2) );
 482 }
 483
 484 // //=======================================================================
 485 // //function : countShape
 486 // //purpose  :
 487 // //=======================================================================
 488 //
 489 // template < class Mesh, class Shape >
 490 // static int countShape( Mesh* mesh, Shape shape ) {
 491 //   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
 492 //   TopTools_MapOfShape mapShape;
 493 //   int nbShape = 0;
 494 //   for ( ; expShape.More(); expShape.Next() ) {
 495 //     if (mapShape.Add(expShape.Current())) {
 496 //       nbShape++;
 497 //     }
 498 //   }
 499 //   return nbShape;
 500 // }
 501 //
 502 // //=======================================================================
 503 // //function : getShape
 504 // //purpose  :
 505 // //=======================================================================
 506 //
 507 // template < class Mesh, class Shape, class Tab >
 508 // void getShape(Mesh* mesh, Shape shape, Tab *t_Shape) {
 509 //   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
 510 //   TopTools_MapOfShape mapShape;
 511 //   for ( int i=0; expShape.More(); expShape.Next() ) {
 512 //     if (mapShape.Add(expShape.Current())) {
 513 //       t_Shape[i] = expShape.Current();
 514 //       i++;
 515 //     }
 516 //   }
 517 //   return;
 518 // }
 519 //
 520 // // //=======================================================================
 521 // // //function : findEdgeID
 522 // // //purpose  :
 523 // // //=======================================================================
 524 //
 525 // static int findEdgeID(const SMDS_MeshNode* aNode,
 526 //                       const SMESHDS_Mesh*  theMesh,
 527 //                       const int            nEdge,
 528 //                       const TopoDS_Shape*  t_Edge) {
 529 //
 530 //   TopoDS_Shape aPntShape, foundEdge;
 531 //   TopoDS_Vertex aVertex;
 532 //   gp_Pnt aPnt( aNode->X(), aNode->Y(), aNode->Z() );
 533 //
 534 //   int foundInd, ind;
 535 //   double nearest = RealLast(), *t_Dist;
 536 //   double epsilon = Precision::Confusion();
 537 //
 538 //   t_Dist = new double[ nEdge ];
 539 //   aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape();
 540 //   aVertex   = TopoDS::Vertex( aPntShape );
 541 //
 542 //   for ( ind=0; ind < nEdge; ind++ ) {
 543 //     BRepExtrema_DistShapeShape aDistance ( aVertex, t_Edge[ind] );
 544 //     t_Dist[ind] = aDistance.Value();
 545 //     if ( t_Dist[ind] < nearest ) {
 546 //       nearest   = t_Dist[ind];
 547 //       foundEdge = t_Edge[ind];
 548 //       foundInd  = ind;
 549 //       if ( nearest < epsilon )
 550 //         ind = nEdge;
 551 //     }
 552 //   }
 553 //
 554 //   delete [] t_Dist;
 555 //   return theMesh->ShapeToIndex( foundEdge );
 556 // }
 557 //
 558 //
 559 // // =======================================================================
 560 // // function : readGMFFile
 561 // // purpose  : read GMF file with geometry associated to mesh
 562 // // =======================================================================
 563 //
 564 // static bool readGMFFile(const int                       fileOpen,
 565 //                         const char*                     theFileName,
 566 //                         SMESH_Mesh&                     theMesh,
 567 //                         const int                       nbShape,
 568 //                         const TopoDS_Shape*             tabShape,
 569 //                         double**                        tabBox,
 570 //                         map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
 571 //                         bool                            toMeshHoles,
 572 //                         int                             nbEnforcedVertices,
 573 //                         int                             nbEnforcedNodes)
 574 // {
 575 //   TopoDS_Shape aShape;
 576 //   TopoDS_Vertex aVertex;
 577 //   SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
 578 //   int nbElem = 0, nbRef = 0, IdShapeRef = 1;
 579 //   int *tabID;
 580 //   int aGMFNodeID = 0;
 581 //   int compoundID =
 582 //     nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
 583 //   int tetraShapeID = compoundID;
 584 //   double epsilon = Precision::Confusion();
 585 //   int *nodeAssigne, *GMFNodeAssigne;
 586 //   SMDS_MeshNode** GMFNode;
 587 //   TopoDS_Shape *tabCorner, *tabEdge;
 588 //   std::map <GmfKwdCod,int> tabRef;
 589 //
 590 //
 591 //   int ver, dim;
 592 //   MESSAGE("Read " << theFileName << " file");
 593 //   int InpMsh = GmfOpenMesh(theFileName, GmfRead, &ver, &dim);
 594 //   if (!InpMsh)
 595 //     return false;
 596 //
 597 //   // ===========================
 598 //   // Fill the tabID array: BEGIN
 599 //   // ===========================
 600 //
 601 //   /*
 602 //   The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
 603 //   */
 604 //   Kernel_Utils::Localizer loc;
 605 //   struct stat status;
 606 //   size_t      length;
 607 //
 608 //   char *ptr, *mapPtr;
 609 //   char *tetraPtr;
 610 //   int *tab = new int[3];
 611 //
 612 //   // Read the file state
 613 //   fstat(fileOpen, &status);
 614 //   length   = status.st_size;
 615 //
 616 //   // Mapping the result file into memory
 617 // #ifdef WNT
 618 //   HANDLE fd = CreateFile(theFileName, GENERIC_READ, FILE_SHARE_READ,
 619 //                          NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
 620 //   HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
 621 //                                         0, (DWORD)length, NULL);
 622 //   ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
 623 // #else
 624 //   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
 625 // #endif
 626 //   mapPtr = ptr;
 627 //
 628 //   ptr      = readMapIntLine(ptr, tab);
 629 //   tetraPtr = ptr;
 630 //
 631 //   nbElem            = tab[0];
 632 //   int nbNodes       = tab[1];
 633 //
 634 //   for (int i=0; i < 4*nbElem; i++)
 635 //     strtol(ptr, &ptr, 10);
 636 //
 637 //   for (int iNode=1; iNode <= nbNodes; iNode++)
 638 //     for (int iCoor=0; iCoor < 3; iCoor++)
 639 //       strtod(ptr, &ptr);
 640 //
 641 //
 642 //   // Reading the number of triangles which corresponds to the number of sub-domains
 643 //   int nbTriangle = strtol(ptr, &ptr, 10);
 644 //
 645 //
 646 //   // The keyword does not exist yet => to update when it is created
 647 // //   int nbSubdomains = GmfStatKwd(InpMsh, GmfSubdomain);
 648 // //   int id_tri[3];
 649 //
 650 //
 651 //   tabID = new int[nbTriangle];
 652 //   for (int i=0; i < nbTriangle; i++) {
 653 //     tabID[i] = 0;
 654 //     int nodeId1, nodeId2, nodeId3;
 655 //     // find the solid corresponding to GHS3D sub-domain following
 656 //     // the technique proposed in GHS3D manual in chapter
 657 //     // "B.4 Subdomain (sub-region) assignment"
 658 //
 659 //     nodeId1 = strtol(ptr, &ptr, 10);
 660 //     nodeId2 = strtol(ptr, &ptr, 10);
 661 //     nodeId3 = strtol(ptr, &ptr, 10);
 662 //
 663 // //   // The keyword does not exist yet => to update when it is created
 664 // //     GmfGetLin(InpMsh, GmfSubdomain, &id_tri[0], &id_tri[1], &id_tri[2]);
 665 // //     nodeId1 = id_tri[0];
 666 // //     nodeId2 = id_tri[1];
 667 // //     nodeId3 = id_tri[2];
 668 //
 669 //     if ( nbTriangle > 1 ) {
 670 //       // get the nodes indices
 671 //       const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
 672 //       const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
 673 //       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
 674 //       try {
 675 //         OCC_CATCH_SIGNALS;
 676 //         tabID[i] = findShapeID( theMesh, n1, n2, n3, toMeshHoles );
 677 //         // -- 0020330: Pb with ghs3d as a submesh
 678 //         // check that found shape is to be meshed
 679 //         if ( tabID[i] > 0 ) {
 680 //           const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] );
 681 //           bool isToBeMeshed = false;
 682 //           for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS )
 683 //             isToBeMeshed = foundShape.IsSame( tabShape[ iS ]);
 684 //           if ( !isToBeMeshed )
 685 //             tabID[i] = HOLE_ID;
 686 //         }
 687 //         // END -- 0020330: Pb with ghs3d as a submesh
 688 // #ifdef _DEBUG_
 689 //         std::cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << std::endl;
 690 // #endif
 691 //       }
 692 //       catch ( Standard_Failure & ex)
 693 //       {
 694 // #ifdef _DEBUG_
 695 //         std::cout << i+1 << " subdomain: Exception caugt: " << ex.GetMessageString() << std::endl;
 696 // #endif
 697 //       }
 698 //       catch (...) {
 699 // #ifdef _DEBUG_
 700 //         std::cout << i+1 << " subdomain: unknown exception caught " << std::endl;
 701 // #endif
 702 //       }
 703 //     }
 704 //   }
 705 //
 706 //   // ===========================
 707 //   // Fill the tabID array: END
 708 //   // ===========================
 709 //
 710 //
 711 //   tabRef[GmfVertices]       = 3;
 712 //   tabRef[GmfCorners]        = 1;
 713 //   tabRef[GmfEdges]          = 2;
 714 //   tabRef[GmfRidges]         = 1;
 715 //   tabRef[GmfTriangles]      = 3;
 716 // //   tabRef[GmfQuadrilaterals] = 4;
 717 //   tabRef[GmfTetrahedra]     = 4;
 718 // //   tabRef[GmfHexahedra]      = 8;
 719 //
 720 //   SMDS_NodeIteratorPtr itOnGMFInputNode = theMeshDS->nodesIterator();
 721 //   while ( itOnGMFInputNode->more() )
 722 //     theMeshDS->RemoveNode( itOnGMFInputNode->next() );
 723 //
 724 //
 725 //   int nbVertices = GmfStatKwd(InpMsh, GmfVertices);
 726 //   int nbCorners = max(countShape( theMeshDS, TopAbs_VERTEX ) , GmfStatKwd(InpMsh, GmfCorners));
 727 //   int nbShapeEdge = countShape( theMeshDS, TopAbs_EDGE );
 728 //
 729 //   tabCorner       = new TopoDS_Shape[ nbCorners ];
 730 //   tabEdge         = new TopoDS_Shape[ nbShapeEdge ];
 731 //   nodeAssigne     = new int[ nbVertices + 1 ];
 732 //   GMFNodeAssigne  = new int[ nbVertices + 1 ];
 733 //   GMFNode         = new SMDS_MeshNode*[ nbVertices + 1 ];
 734 //
 735 //   getShape(theMeshDS, TopAbs_VERTEX, tabCorner);
 736 //   getShape(theMeshDS, TopAbs_EDGE,   tabEdge);
 737 //
 738 //   std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
 739 //   for ( ; it != tabRef.end() ; ++it)
 740 //   {
 741 // //     int dummy;
 742 //     GmfKwdCod token = it->first;
 743 //     nbRef    = it->second;
 744 //
 745 //     nbElem = GmfStatKwd(InpMsh, token);
 746 //     if (nbElem > 0) {
 747 //       GmfGotoKwd(InpMsh, token);
 748 //       std::cout << "Read " << nbElem;
 749 //     }
 750 //     else
 751 //       continue;
 752 //
 753 //     int id[nbElem*tabRef[token]];
 754 //     int ghs3dShapeID[nbElem];
 755 //
 756 //     if (token == GmfVertices) {
 757 //       std::cout << " vertices" << std::endl;
 758 //       int aGMFID;
 759 //
 760 //       float VerTab_f[nbElem][3];
 761 //       double VerTab_d[nbElem][3];
 762 //       SMDS_MeshNode * aGMFNode;
 763 //
 764 //       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
 765 //         aGMFID = iElem + 1;
 766 //         if (ver == GmfFloat) {
 767 //           GmfGetLin(InpMsh, token, &VerTab_f[nbElem][0], &VerTab_f[nbElem][1], &VerTab_f[nbElem][2], &ghs3dShapeID[iElem]);
 768 //           aGMFNode = theMeshDS->AddNode(VerTab_f[nbElem][0], VerTab_f[nbElem][1], VerTab_f[nbElem][2]);
 769 //         }
 770 //         else {
 771 //           GmfGetLin(InpMsh, token, &VerTab_d[nbElem][0], &VerTab_d[nbElem][1], &VerTab_d[nbElem][2], &ghs3dShapeID[iElem]);
 772 //           aGMFNode = theMeshDS->AddNode(VerTab_d[nbElem][0], VerTab_d[nbElem][1], VerTab_d[nbElem][2]);
 773 //         }
 774 //         GMFNode[ aGMFID ] = aGMFNode;
 775 //         nodeAssigne[ aGMFID ] = 0;
 776 //         GMFNodeAssigne[ aGMFID ] = 0;
 777 //       }
 778 //     }
 779 //     else if (token == GmfCorners && nbElem > 0) {
 780 //       std::cout << " corners" << std::endl;
 781 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 782 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
 783 //     }
 784 //     else if (token == GmfRidges && nbElem > 0) {
 785 //       std::cout << " ridges" << std::endl;
 786 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 787 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
 788 //     }
 789 //     else if (token == GmfEdges && nbElem > 0) {
 790 //       std::cout << " edges" << std::endl;
 791 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 792 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &ghs3dShapeID[iElem]);
 793 //     }
 794 //     else if (token == GmfTriangles && nbElem > 0) {
 795 //       std::cout << " triangles" << std::endl;
 796 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 797 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &ghs3dShapeID[iElem]);
 798 //     }
 799 // //     else if (token == GmfQuadrilaterals && nbElem > 0) {
 800 // //       std::cout << " Quadrilaterals" << std::endl;
 801 // //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 802 // //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &ghs3dShapeID[iElem]);
 803 // //     }
 804 //     else if (token == GmfTetrahedra && nbElem > 0) {
 805 //       std::cout << " Tetrahedra" << std::endl;
 806 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 807 //         GmfGetLin(InpMsh, token,
 808 //                   &id[iElem*tabRef[token]],
 809 //                   &id[iElem*tabRef[token]+1],
 810 //                   &id[iElem*tabRef[token]+2],
 811 //                   &id[iElem*tabRef[token]+3],
 812 //                   &ghs3dShapeID[iElem]);
 813 //     }
 814 // //     else if (token == GmfHexahedra && nbElem > 0) {
 815 // //       std::cout << " Hexahedra" << std::endl;
 816 // //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 817 // //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
 818 // //                   &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &ghs3dShapeID[iElem]);
 819 // //     }
 820 //
 821 //     switch (token) {
 822 //     case GmfCorners:
 823 //     case GmfRidges:
 824 //     case GmfEdges:
 825 //     case GmfTriangles:
 826 // //     case GmfQuadrilaterals:
 827 //     case GmfTetrahedra:
 828 // //     case GmfHexahedra:
 829 //     {
 830 //       int nodeDim, shapeID, *nodeID;
 831 //       const SMDS_MeshNode** node;
 832 // //       std::vector< SMDS_MeshNode* > enfNode( nbRef );
 833 //       SMDS_MeshElement * aGMFElement;
 834 //
 835 //       node    = new const SMDS_MeshNode*[nbRef];
 836 //       nodeID  = new int[ nbRef ];
 837 //
 838 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 839 //       {
 840 //         for ( int iRef = 0; iRef < nbRef; iRef++ )
 841 //         {
 842 //           aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
 843 //           node  [ iRef ] = GMFNode[ aGMFNodeID ];
 844 //           nodeID[ iRef ] = aGMFNodeID;
 845 //         }
 846 //
 847 //         switch (token)
 848 //         {
 849 //         case GmfCorners: {
 850 //           nodeDim = 1;
 851 //           gp_Pnt GMFPnt ( node[0]->X(), node[0]->Y(), node[0]->Z() );
 852 //           for ( int i=0; i<nbElem; i++ ) {
 853 //             aVertex = TopoDS::Vertex( tabCorner[i] );
 854 //             gp_Pnt aPnt = BRep_Tool::Pnt( aVertex );
 855 //             if ( aPnt.Distance( GMFPnt ) < epsilon )
 856 //               break;
 857 //           }
 858 //           break;
 859 //         }
 860 //         case GmfEdges: {
 861 //           nodeDim = 2;
 862 //           aGMFElement = theMeshDS->AddEdge( node[0], node[1] );
 863 //           int iNode = 1;
 864 //           if ( GMFNodeAssigne[ nodeID[0] ] == 0 || GMFNodeAssigne[ nodeID[0] ] == 2 )
 865 //             iNode = 0;
 866 //           shapeID = findEdgeID( node[iNode], theMeshDS, nbShapeEdge, tabEdge );
 867 //           break;
 868 //         }
 869 //         case GmfRidges:
 870 //           break;
 871 //         case GmfTriangles: {
 872 //           nodeDim = 3;
 873 //           aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2]);
 874 //           shapeID = -1;
 875 //           break;
 876 //         }
 877 // //         case GmfQuadrilaterals: {
 878 // //           nodeDim = 4;
 879 // //           aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2], node[3] );
 880 // //           shapeID = -1;
 881 // //           break;
 882 // //         }
 883 //         case GmfTetrahedra: {
 884 //
 885 //           // IN WORK
 886 //           TopoDS_Shape aSolid;
 887 //           // We always run GHS3D with "to mesh holes"==TRUE but we must not create
 888 //           // tetras within holes depending on hypo option,
 889 //           // so we first check if aTet is inside a hole and then create it
 890 //           if ( nbTriangle > 1 ) {
 891 //             tetraShapeID = HOLE_ID; // negative tetraShapeID means not to create tetras if !toMeshHoles
 892 //             int aGhs3dShapeID = ghs3dShapeID[iElem] - IdShapeRef;
 893 //             if ( tabID[ aGhs3dShapeID ] == 0 ) {
 894 //               TopAbs_State state;
 895 //               aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
 896 //               if ( toMeshHoles || state == TopAbs_IN )
 897 //                 tetraShapeID = theMeshDS->ShapeToIndex( aSolid );
 898 //               tabID[ aGhs3dShapeID ] = tetraShapeID;
 899 //             }
 900 //             else
 901 //               tetraShapeID = tabID[ aGhs3dShapeID ];
 902 //           }
 903 //           else if ( nbShape > 1 ) {
 904 //             // Case where nbTriangle == 1 while nbShape == 2 encountered
 905 //             // with compound of 2 boxes and "To mesh holes"==False,
 906 //             // so there are no subdomains specified for each tetrahedron.
 907 //             // Try to guess a solid by a node already bound to shape
 908 //             tetraShapeID = 0;
 909 //             for ( int i=0; i<4 && tetraShapeID==0; i++ ) {
 910 //               if ( nodeAssigne[ nodeID[i] ] == 1 &&
 911 //                   node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
 912 //                   node[i]->getshapeId() > 1 )
 913 //               {
 914 //                 tetraShapeID = node[i]->getshapeId();
 915 //               }
 916 //             }
 917 //             if ( tetraShapeID==0 ) {
 918 //               aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
 919 //               tetraShapeID = theMeshDS->ShapeToIndex( aSolid );
 920 //             }
 921 //           }
 922 //           // set new nodes and tetrahedron onto the shape
 923 //           for ( int i=0; i<4; i++ ) {
 924 //             if ( nodeAssigne[ nodeID[i] ] == 0 ) {
 925 //               if ( tetraShapeID != HOLE_ID )
 926 //                 theMeshDS->SetNodeInVolume( node[i], tetraShapeID );
 927 //               nodeAssigne[ nodeID[i] ] = tetraShapeID;
 928 //             }
 929 //           }
 930 //           if ( toMeshHoles || tetraShapeID != HOLE_ID ) {
 931 //             aGMFElement = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
 932 //             theMeshDS->SetMeshElementOnShape( aGMFElement, tetraShapeID );
 933 //           }
 934 //
 935 //           // IN WORK
 936 //
 937 //           nodeDim = 5;
 938 //           break;
 939 //         }
 940 // //         case GmfHexahedra: {
 941 // //           nodeDim = 6;
 942 // //           aGMFElement = theMeshDS->AddVolume( node[0], node[3], node[2], node[1],
 943 // //                                             node[4], node[7], node[6], node[5] );
 944 // //           break;
 945 // //         }
 946 //         default: continue;
 947 //         }
 948 //         if (token != GmfRidges)
 949 //         {
 950 //           for ( int i=0; i<nbRef; i++ ) {
 951 //               if ( GMFNodeAssigne[ nodeID[i] ] == 0 ) {
 952 //                 if      ( token == GmfCorners )   theMeshDS->SetNodeOnVertex( node[0], aVertex );
 953 //                 else if ( token == GmfEdges )     theMeshDS->SetNodeOnEdge( node[i], shapeID );
 954 //                 else if ( token == GmfTriangles ) theMeshDS->SetNodeOnFace( node[i], shapeID );
 955 //                 GMFNodeAssigne[ nodeID[i] ] = nodeDim;
 956 //               }
 957 //             }
 958 //             if ( token != "Corners" )
 959 //               theMeshDS->SetMeshElementOnShape( aGMFElement, shapeID );
 960 //         }
 961 //       } // for
 962 //
 963 //       if ( !toMeshHoles ) {
 964 //         map <int,const SMDS_MeshNode*>::iterator itOnNode = theGhs3dIdToNodeMap.find( nbVertices-(nbEnforcedVertices+nbEnforcedNodes) );
 965 //         for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
 966 //           if ( nodeAssigne[ itOnNode->first ] == HOLE_ID )
 967 //             theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
 968 //         }
 969 //       }
 970 //
 971 //       delete [] node;
 972 //       delete [] nodeID;
 973 //       break;
 974 //       } // case GmfTetrahedra
 975 //     } // switch(token)
 976 //   } // for
 977 //   cout << std::endl;
 978 //
 979 // #ifdef WNT
 980 //   UnmapViewOfFile(mapPtr);
 981 //   CloseHandle(hMapObject);
 982 //   CloseHandle(fd);
 983 // #else
 984 //   munmap(mapPtr, length);
 985 // #endif
 986 //   close(fileOpen);
 987 //
 988 //   delete [] tabID;
 989 //   delete [] tabCorner;
 990 //   delete [] tabEdge;
 991 //   delete [] nodeAssigne;
 992 //   delete [] GMFNodeAssigne;
 993 //   delete [] GMFNode;
 994 //
 995 //   return true;
 996 // }
 997
 998
 999 //=======================================================================
1000 //function : addElemInMeshGroup
1001 //purpose  : Update or create groups in mesh
1002 //=======================================================================
1003
1004 static void addElemInMeshGroup(SMESH_Mesh*             theMesh,
1005                                const SMDS_MeshElement* anElem,
1006                                std::string&            groupName,
1007                                std::set<std::string>&  groupsToRemove)
1008 {
1009   if ( !anElem ) return; // issue 0021776
1010
1011   bool groupDone = false;
1012   SMESH_Mesh::GroupIteratorPtr grIt = theMesh->GetGroups();
1013   while (grIt->more()) {
1014     SMESH_Group * group = grIt->next();
1015     if ( !group ) continue;
1016     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
1017     if ( !groupDS ) continue;
1018     if ( groupDS->GetType()==anElem->GetType() &&groupName.compare(group->GetName())==0) {
1019       SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
1020       aGroupDS->SMDSGroup().Add(anElem);
1021       groupDone = true;
1022 //       MESSAGE("Successfully added enforced element to existing group " << groupName);
1023       break;
1024     }
1025   }
1026
1027   if (!groupDone)
1028   {
1029     int groupId;
1030     SMESH_Group* aGroup = theMesh->AddGroup(anElem->GetType(), groupName.c_str(), groupId);
1031     aGroup->SetName( groupName.c_str() );
1032     SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
1033     aGroupDS->SMDSGroup().Add(anElem);
1034 //     MESSAGE("Successfully created enforced vertex group " << groupName);
1035     groupDone = true;
1036   }
1037   if (!groupDone)
1038     throw SALOME_Exception(LOCALIZED("A given element was not added to a group"));
1039 }
1040
1041
1042 //=======================================================================
1043 //function : updateMeshGroups
1044 //purpose  : Update or create groups in mesh
1045 //=======================================================================
1046
1047 static void updateMeshGroups(SMESH_Mesh* theMesh, std::set<std::string> groupsToRemove)
1048 {
1049   SMESH_Mesh::GroupIteratorPtr grIt = theMesh->GetGroups();
1050   while (grIt->more()) {
1051     SMESH_Group * group = grIt->next();
1052     if ( !group ) continue;
1053     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
1054     if ( !groupDS ) continue;
1055     std::string currentGroupName = (string)group->GetName();
1056     if (groupDS->IsEmpty() && groupsToRemove.find(currentGroupName) != groupsToRemove.end()) {
1057       // Previous group created by enforced elements
1058       MESSAGE("Delete previous group created by removed enforced elements: " << group->GetName())
1059       theMesh->RemoveGroup(groupDS->GetID());
1060     }
1061   }
1062 }
1063
1064 //=======================================================================
1065 //function : removeEmptyGroupsOfDomains
1066 //purpose  : remove empty groups named "Domain_nb" created due to
1067 //           "To make groups of domains" option.
1068 //=======================================================================
1069
1070 static void removeEmptyGroupsOfDomains(SMESH_Mesh* mesh,
1071                                        bool        notEmptyAsWell = false)
1072 {
1073   const char* refName = theDomainGroupNamePrefix;
1074   const size_t refLen = strlen( theDomainGroupNamePrefix );
1075
1076   std::list<int> groupIDs = mesh->GetGroupIds();
1077   std::list<int>::const_iterator id = groupIDs.begin();
1078   for ( ; id != groupIDs.end(); ++id )
1079   {
1080     SMESH_Group* group = mesh->GetGroup( *id );
1081     if ( !group || ( !group->GetGroupDS()->IsEmpty() && !notEmptyAsWell ))
1082       continue;
1083     const char* name = group->GetName();
1084     char* end;
1085     // check the name
1086     if ( strncmp( name, refName, refLen ) == 0 && // starts from refName;
1087          isdigit( *( name + refLen )) &&          // refName is followed by a digit;
1088          strtol( name + refLen, &end, 10) >= 0 && // there are only digits ...
1089          *end == '\0')                            // ... till a string end.
1090     {
1091       mesh->RemoveGroup( *id );
1092     }
1093   }
1094 }
1095
1096 //================================================================================
1097 /*!
1098  * \brief Create the groups corresponding to domains
1099  */
1100 //================================================================================
1101
1102 static void makeDomainGroups( std::vector< std::vector< const SMDS_MeshElement* > >& elemsOfDomain,
1103                               SMESH_MesherHelper*                                    theHelper)
1104 {
1105   // int nbDomains = 0;
1106   // for ( size_t i = 0; i < elemsOfDomain.size(); ++i )
1107   //   nbDomains += ( elemsOfDomain[i].size() > 0 );
1108
1109   // if ( nbDomains > 1 )
1110   for ( size_t iDomain = 0; iDomain < elemsOfDomain.size(); ++iDomain )
1111   {
1112     std::vector< const SMDS_MeshElement* > & elems = elemsOfDomain[ iDomain ];
1113     if ( elems.empty() ) continue;
1114
1115     // find existing groups
1116     std::vector< SMESH_Group* > groupOfType( SMDSAbs_NbElementTypes, (SMESH_Group*)NULL );
1117     const std::string domainName = ( SMESH_Comment( theDomainGroupNamePrefix ) << iDomain );
1118     SMESH_Mesh::GroupIteratorPtr groupIt = theHelper->GetMesh()->GetGroups();
1119     while ( groupIt->more() )
1120     {
1121       SMESH_Group* group = groupIt->next();
1122       if ( domainName == group->GetName() &&
1123            dynamic_cast< SMESHDS_Group* >( group->GetGroupDS()) )
1124         groupOfType[ group->GetGroupDS()->GetType() ] = group;
1125     }
1126     // create and fill the groups
1127     size_t iElem = 0;
1128     int groupID;
1129     do
1130     {
1131       SMESH_Group* group = groupOfType[ elems[ iElem ]->GetType() ];
1132       if ( !group )
1133         group = theHelper->GetMesh()->AddGroup( elems[ iElem ]->GetType(),
1134                                                 domainName.c_str(), groupID );
1135       SMDS_MeshGroup& groupDS =
1136         static_cast< SMESHDS_Group* >( group->GetGroupDS() )->SMDSGroup();
1137
1138       while ( iElem < elems.size() && groupDS.Add( elems[iElem] ))
1139         ++iElem;
1140
1141     } while ( iElem < elems.size() );
1142   }
1143 }
1144
1145 //=======================================================================
1146 //function : readGMFFile
1147 //purpose  : read GMF file w/o geometry associated to mesh
1148 //=======================================================================
1149
1150 static bool readGMFFile(const char*                     theFile,
1151                         GHS3DPlugin_GHS3D*              theAlgo,
1152                         SMESH_MesherHelper*             theHelper,
1153                         TopoDS_Shape                    theSolid,
1154                         vector <const SMDS_MeshNode*> & theNodeByGhs3dId,
1155                         map<const SMDS_MeshNode*,int> & theNodeToGhs3dIdMap,
1156                         std::vector<std::string> &      aNodeGroupByGhs3dId,
1157                         std::vector<std::string> &      anEdgeGroupByGhs3dId,
1158                         std::vector<std::string> &      aFaceGroupByGhs3dId,
1159                         std::set<std::string> &         groupsToRemove,
1160                         bool                            toMakeGroupsOfDomains=false)
1161 {
1162   std::string tmpStr;
1163   SMESHDS_Mesh* theMeshDS = theHelper->GetMeshDS();
1164
1165   int nbInitialNodes = theNodeByGhs3dId.size();
1166   int nbMeshNodes = theMeshDS->NbNodes();
1167
1168   const bool isQuadMesh =
1169     theHelper->GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
1170     theHelper->GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
1171     theHelper->GetMesh()->NbVolumes( ORDER_QUADRATIC );
1172
1173 #ifdef _DEBUG_
1174   std::cout << "theNodeByGhs3dId.size(): " << nbInitialNodes << std::endl;
1175   std::cout << "theHelper->GetMesh()->NbNodes(): " << nbMeshNodes << std::endl;
1176   std::cout << "isQuadMesh: " << isQuadMesh << std::endl;
1177 #endif
1178
1179   if (theHelper->GetSubShapeID() != 0)
1180     theHelper->IsQuadraticSubMesh( theHelper->GetSubShape() );
1181
1182   // ---------------------------------
1183   // Read generated elements and nodes
1184   // ---------------------------------
1185
1186   int nbElem = 0, nbRef = 0;
1187   int aGMFNodeID = 0;
1188   const SMDS_MeshNode** GMFNode;
1189 #ifdef _DEBUG_
1190   std::map<int, std::set<int> > subdomainId2tetraId;
1191 #endif
1192   std::map <GmfKwdCod,int> tabRef;
1193   const bool force3d = true; // since there is no geometry
1194   const int  noID  = 0;
1195
1196   tabRef[GmfVertices]       = 3; // for new nodes and enforced nodes
1197   tabRef[GmfCorners]        = 1;
1198   tabRef[GmfEdges]          = 2; // for enforced edges
1199   tabRef[GmfRidges]         = 1;
1200   tabRef[GmfTriangles]      = 3; // for enforced faces
1201   tabRef[GmfQuadrilaterals] = 4;
1202   tabRef[GmfTetrahedra]     = 4; // for new tetras
1203   tabRef[GmfHexahedra]      = 8;
1204
1205   int ver, dim;
1206   MESSAGE("Read " << theFile << " file");
1207   int InpMsh = GmfOpenMesh(theFile, GmfRead, &ver, &dim);
1208   if (!InpMsh)
1209     return false;
1210   MESSAGE("Done ");
1211
1212   // Issue 0020682. Avoid creating nodes and tetras at place where
1213   // volumic elements already exist
1214   SMESH_ElementSearcher* elemSearcher = 0;
1215   std::vector< const SMDS_MeshElement* > foundVolumes;
1216   if ( theHelper->GetMesh()->NbVolumes() > 0 )
1217     elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *theHelper->GetMeshDS() );
1218
1219   // IMP 0022172: [CEA 790] create the groups corresponding to domains
1220   std::vector< std::vector< const SMDS_MeshElement* > > elemsOfDomain;
1221
1222   int nbVertices = GmfStatKwd(InpMsh, GmfVertices) - nbInitialNodes;
1223   GMFNode = new const SMDS_MeshNode*[ nbVertices + 1 ];
1224
1225   std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
1226   for ( ; it != tabRef.end() ; ++it)
1227   {
1228     if(theAlgo->computeCanceled()) {
1229       GmfCloseMesh(InpMsh);
1230       delete [] GMFNode;
1231       return false;
1232     }
1233     int dummy;
1234     GmfKwdCod token = it->first;
1235     nbRef           = it->second;
1236
1237     nbElem = GmfStatKwd(InpMsh, token);
1238     if (nbElem > 0) {
1239       GmfGotoKwd(InpMsh, token);
1240       std::cout << "Read " << nbElem;
1241     }
1242     else
1243       continue;
1244
1245     std::vector<int> id (nbElem*tabRef[token]); // node ids
1246     std::vector<int> domainID( nbElem ); // domain
1247
1248     if (token == GmfVertices) {
1249       (nbElem <= 1) ? tmpStr = " vertex" : tmpStr = " vertices";
1250 //       std::cout << nbInitialNodes << " from input mesh " << std::endl;
1251
1252       // Remove orphan nodes from previous enforced mesh which was cleared
1253 //       if ( nbElem < nbMeshNodes ) {
1254 //         const SMDS_MeshNode* node;
1255 //         SMDS_NodeIteratorPtr nodeIt = theMeshDS->nodesIterator();
1256 //         while ( nodeIt->more() )
1257 //         {
1258 //           node = nodeIt->next();
1259 //           if (theNodeToGhs3dIdMap.find(node) != theNodeToGhs3dIdMap.end())
1260 //             theMeshDS->RemoveNode(node);
1261 //         }
1262 //       }
1263
1264
1265       int aGMFID;
1266
1267       float VerTab_f[3];
1268       double x, y, z;
1269       const SMDS_MeshNode * aGMFNode;
1270
1271       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
1272         if(theAlgo->computeCanceled()) {
1273           GmfCloseMesh(InpMsh);
1274           delete [] GMFNode;
1275           return false;
1276         }
1277         if (ver == GmfFloat) {
1278           GmfGetLin(InpMsh, token, &VerTab_f[0], &VerTab_f[1], &VerTab_f[2], &dummy);
1279           x = VerTab_f[0];
1280           y = VerTab_f[1];
1281           z = VerTab_f[2];
1282         }
1283         else {
1284           GmfGetLin(InpMsh, token, &x, &y, &z, &dummy);
1285         }
1286         if (iElem >= nbInitialNodes) {
1287           if ( elemSearcher &&
1288                 elemSearcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_Volume, foundVolumes))
1289             aGMFNode = 0;
1290           else
1291             aGMFNode = theHelper->AddNode(x, y, z);
1292
1293           aGMFID = iElem -nbInitialNodes +1;
1294           GMFNode[ aGMFID ] = aGMFNode;
1295           if (aGMFID-1 < aNodeGroupByGhs3dId.size() && !aNodeGroupByGhs3dId.at(aGMFID-1).empty())
1296             addElemInMeshGroup(theHelper->GetMesh(), aGMFNode, aNodeGroupByGhs3dId.at(aGMFID-1), groupsToRemove);
1297         }
1298       }
1299     }
1300     else if (token == GmfCorners && nbElem > 0) {
1301       (nbElem <= 1) ? tmpStr = " corner" : tmpStr = " corners";
1302       for ( int iElem = 0; iElem < nbElem; iElem++ )
1303         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
1304     }
1305     else if (token == GmfRidges && nbElem > 0) {
1306       (nbElem <= 1) ? tmpStr = " ridge" : tmpStr = " ridges";
1307       for ( int iElem = 0; iElem < nbElem; iElem++ )
1308         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
1309     }
1310     else if (token == GmfEdges && nbElem > 0) {
1311       (nbElem <= 1) ? tmpStr = " edge" : tmpStr = " edges";
1312       for ( int iElem = 0; iElem < nbElem; iElem++ )
1313         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &domainID[iElem]);
1314     }
1315     else if (token == GmfTriangles && nbElem > 0) {
1316       (nbElem <= 1) ? tmpStr = " triangle" : tmpStr = " triangles";
1317       for ( int iElem = 0; iElem < nbElem; iElem++ )
1318         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &domainID[iElem]);
1319     }
1320     else if (token == GmfQuadrilaterals && nbElem > 0) {
1321       (nbElem <= 1) ? tmpStr = " Quadrilateral" : tmpStr = " Quadrilaterals";
1322       for ( int iElem = 0; iElem < nbElem; iElem++ )
1323         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &domainID[iElem]);
1324     }
1325     else if (token == GmfTetrahedra && nbElem > 0) {
1326       (nbElem <= 1) ? tmpStr = " Tetrahedron" : tmpStr = " Tetrahedra";
1327       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
1328         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &domainID[iElem]);
1329 #ifdef _DEBUG_
1330         subdomainId2tetraId[dummy].insert(iElem+1);
1331 //         MESSAGE("subdomainId2tetraId["<<dummy<<"].insert("<<iElem+1<<")");
1332 #endif
1333       }
1334     }
1335     else if (token == GmfHexahedra && nbElem > 0) {
1336       (nbElem <= 1) ? tmpStr = " Hexahedron" : tmpStr = " Hexahedra";
1337       for ( int iElem = 0; iElem < nbElem; iElem++ )
1338         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
1339                   &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &domainID[iElem]);
1340     }
1341     std::cout << tmpStr << std::endl;
1342     std::cout << std::endl;
1343
1344     switch (token) {
1345     case GmfCorners:
1346     case GmfRidges:
1347     case GmfEdges:
1348     case GmfTriangles:
1349     case GmfQuadrilaterals:
1350     case GmfTetrahedra:
1351     case GmfHexahedra:
1352     {
1353       std::vector< const SMDS_MeshNode* > node( nbRef );
1354       std::vector< int >          nodeID( nbRef );
1355       std::vector< SMDS_MeshNode* > enfNode( nbRef );
1356       const SMDS_MeshElement* aCreatedElem;
1357
1358       for ( int iElem = 0; iElem < nbElem; iElem++ )
1359       {
1360         if(theAlgo->computeCanceled()) {
1361           GmfCloseMesh(InpMsh);
1362           delete [] GMFNode;
1363           return false;
1364         }
1365         // Check if elem is already in input mesh. If yes => skip
1366         bool fullyCreatedElement = false; // if at least one of the nodes was created
1367         for ( int iRef = 0; iRef < nbRef; iRef++ )
1368         {
1369           aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
1370           if (aGMFNodeID <= nbInitialNodes) // input nodes
1371           {
1372             aGMFNodeID--;
1373             node[ iRef ] = theNodeByGhs3dId[aGMFNodeID];
1374           }
1375           else
1376           {
1377             fullyCreatedElement = true;
1378             aGMFNodeID -= nbInitialNodes;
1379             nodeID[ iRef ] = aGMFNodeID ;
1380             node  [ iRef ] = GMFNode[ aGMFNodeID ];
1381           }
1382         }
1383         aCreatedElem = 0;
1384         switch (token)
1385         {
1386         case GmfEdges:
1387           if (fullyCreatedElement) {
1388             aCreatedElem = theHelper->AddEdge( node[0], node[1], noID, force3d );
1389             if (anEdgeGroupByGhs3dId.size() && !anEdgeGroupByGhs3dId[iElem].empty())
1390               addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, anEdgeGroupByGhs3dId[iElem], groupsToRemove);
1391           }
1392           break;
1393         case GmfTriangles:
1394           if (fullyCreatedElement) {
1395             aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], noID, force3d );
1396             if (aFaceGroupByGhs3dId.size() && !aFaceGroupByGhs3dId[iElem].empty())
1397               addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, aFaceGroupByGhs3dId[iElem], groupsToRemove);
1398           }
1399           break;
1400         case GmfQuadrilaterals:
1401           if (fullyCreatedElement) {
1402             aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], node[3], noID, force3d );
1403           }
1404           break;
1405         case GmfTetrahedra:
1406           if ( elemSearcher ) {
1407             // Issue 0020682. Avoid creating nodes and tetras at place where
1408             // volumic elements already exist
1409             if ( !node[1] || !node[0] || !node[2] || !node[3] )
1410               continue;
1411             if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) +
1412                                                     SMESH_TNodeXYZ(node[1]) +
1413                                                     SMESH_TNodeXYZ(node[2]) +
1414                                                     SMESH_TNodeXYZ(node[3]) ) / 4.,
1415                                                     SMDSAbs_Volume, foundVolumes ))
1416             break;
1417           }
1418           aCreatedElem = theHelper->AddVolume( node[1], node[0], node[2], node[3], noID, force3d );
1419           break;
1420         case GmfHexahedra:
1421           if ( elemSearcher ) {
1422             // Issue 0020682. Avoid creating nodes and tetras at place where
1423             // volumic elements already exist
1424             if ( !node[1] || !node[0] || !node[2] || !node[3] || !node[4] || !node[5] || !node[6] || !node[7])
1425               continue;
1426             if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) +
1427                                                     SMESH_TNodeXYZ(node[1]) +
1428                                                     SMESH_TNodeXYZ(node[2]) +
1429                                                     SMESH_TNodeXYZ(node[3]) +
1430                                                     SMESH_TNodeXYZ(node[4]) +
1431                                                     SMESH_TNodeXYZ(node[5]) +
1432                                                     SMESH_TNodeXYZ(node[6]) +
1433                                                     SMESH_TNodeXYZ(node[7])) / 8.,
1434                                                     SMDSAbs_Volume, foundVolumes ))
1435             break;
1436           }
1437           aCreatedElem = theHelper->AddVolume( node[0], node[3], node[2], node[1],
1438                                                node[4], node[7], node[6], node[5], noID, force3d );
1439           break;
1440         default: continue;
1441         }
1442         if ( aCreatedElem && toMakeGroupsOfDomains )
1443         {
1444           if ( domainID[iElem] >= (int) elemsOfDomain.size() )
1445             elemsOfDomain.resize( domainID[iElem] + 1 );
1446           elemsOfDomain[ domainID[iElem] ].push_back( aCreatedElem );
1447         }
1448       } // loop on elements of one type
1449       break;
1450     } // case ...
1451     } // switch (token)
1452   } // loop on tabRef
1453
1454   GmfCloseMesh(InpMsh);
1455   delete [] GMFNode;
1456
1457   // 0022172: [CEA 790] create the groups corresponding to domains
1458   if ( toMakeGroupsOfDomains )
1459     makeDomainGroups( elemsOfDomain, theHelper );
1460
1461 #ifdef _DEBUG_
1462   MESSAGE("Nb subdomains " << subdomainId2tetraId.size());
1463   std::map<int, std::set<int> >::const_iterator subdomainIt = subdomainId2tetraId.begin();
1464   TCollection_AsciiString aSubdomainFileName = theFile;
1465   aSubdomainFileName = aSubdomainFileName + ".subdomain";
1466   ofstream aSubdomainFile  ( aSubdomainFileName.ToCString()  , ios::out);
1467
1468   aSubdomainFile << "Nb subdomains " << subdomainId2tetraId.size() << std::endl;
1469   for(;subdomainIt != subdomainId2tetraId.end() ; ++subdomainIt) {
1470     int subdomainId = subdomainIt->first;
1471     std::set<int> tetraIds = subdomainIt->second;
1472     MESSAGE("Subdomain #"<<subdomainId<<": "<<tetraIds.size()<<" tetrahedrons");
1473     std::set<int>::const_iterator tetraIdsIt = tetraIds.begin();
1474     aSubdomainFile << subdomainId << std::endl;
1475     for(;tetraIdsIt != tetraIds.end() ; ++tetraIdsIt) {
1476       aSubdomainFile << (*tetraIdsIt) << " ";
1477     }
1478     aSubdomainFile << std::endl;
1479   }
1480   aSubdomainFile.close();
1481 #endif
1482
1483   if (elemSearcher)
1484   {
1485     MESSAGE("delete elemSearcher")
1486     delete elemSearcher;
1487     elemSearcher = 0;
1488   }
1489   return true;
1490 }
1491
1492 static bool writeGMFFile(const char*                                     theMeshFileName,
1493                          const char*                                     theRequiredFileName,
1494                          const char*                                     theSolFileName,
1495                          const SMESH_ProxyMesh&                          theProxyMesh,
1496                          SMESH_Mesh *                                    theMesh,
1497                          std::vector <const SMDS_MeshNode*> &            theNodeByGhs3dId,
1498                          std::map<const SMDS_MeshNode*,int> &            aNodeToGhs3dIdMap,
1499                          std::vector<std::string> &                      aNodeGroupByGhs3dId,
1500                          std::vector<std::string> &                      anEdgeGroupByGhs3dId,
1501                          std::vector<std::string> &                      aFaceGroupByGhs3dId,
1502                          GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap & theEnforcedNodes,
1503                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges,
1504                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles,
1505                          std::map<std::vector<double>, std::string> &    enfVerticesWithGroup,
1506                          GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices)
1507 {
1508   MESSAGE("writeGMFFile w/o geometry");
1509   std::string tmpStr;
1510   int idx, idxRequired = 0, idxSol = 0;
1511   const int dummyint = 0;
1512   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt;
1513   std::vector<double> enfVertexSizes;
1514   const SMDS_MeshElement* elem;
1515   TIDSortedElemSet anElemSet, theKeptEnforcedEdges, theKeptEnforcedTriangles;
1516   SMDS_ElemIteratorPtr nodeIt;
1517   std::vector <const SMDS_MeshNode*> theEnforcedNodeByGhs3dId;
1518   map<const SMDS_MeshNode*,int> anEnforcedNodeToGhs3dIdMap, anExistingEnforcedNodeToGhs3dIdMap;
1519   std::vector< const SMDS_MeshElement* > foundElems;
1520   map<const SMDS_MeshNode*,TopAbs_State> aNodeToTopAbs_StateMap;
1521   int nbFoundElems;
1522   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap::iterator elemIt;
1523   TIDSortedElemSet::iterator elemSetIt;
1524   bool isOK;
1525   auto_ptr< SMESH_ElementSearcher > pntCls
1526     ( SMESH_MeshAlgos::GetElementSearcher(*theMesh->GetMeshDS()));
1527
1528   int nbEnforcedVertices = theEnforcedVertices.size();
1529
1530   // count faces
1531   int nbFaces = theProxyMesh.NbFaces();
1532   int nbNodes;
1533
1534   // groups management
1535   int usedEnforcedNodes = 0;
1536   std::string gn = "";
1537
1538   if ( nbFaces == 0 )
1539     return false;
1540
1541   idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1542   if (!idx)
1543     return false;
1544
1545   /* ========================== FACES ========================== */
1546   /* TRIANGLES ========================== */
1547   SMDS_ElemIteratorPtr eIt = theProxyMesh.GetFaces();
1548   while ( eIt->more() )
1549   {
1550     elem = eIt->next();
1551     anElemSet.insert(elem);
1552     nodeIt = elem->nodesIterator();
1553     nbNodes = elem->NbCornerNodes();
1554     while ( nodeIt->more() && nbNodes--)
1555     {
1556       // find GHS3D ID
1557       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1558       int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
1559       aNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1560     }
1561   }
1562
1563   /* EDGES ========================== */
1564
1565   // Iterate over the enforced edges
1566   for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
1567     elem = elemIt->first;
1568     isOK = true;
1569     nodeIt = elem->nodesIterator();
1570     nbNodes = 2;
1571     while ( nodeIt->more() && nbNodes-- ) {
1572       // find GHS3D ID
1573       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1574       // Test if point is inside shape to mesh
1575       gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1576       TopAbs_State result = pntCls->GetPointState( myPoint );
1577       if ( result == TopAbs_OUT ) {
1578         isOK = false;
1579         break;
1580       }
1581       aNodeToTopAbs_StateMap.insert( make_pair( node, result ));
1582     }
1583     if (isOK) {
1584       nodeIt = elem->nodesIterator();
1585       nbNodes = 2;
1586       int newId = -1;
1587       while ( nodeIt->more() && nbNodes-- ) {
1588         // find GHS3D ID
1589         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1590         gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1591         nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1592 #ifdef _DEBUG_
1593         std::cout << "Node at "<<node->X()<<", "<<node->Y()<<", "<<node->Z()<<std::endl;
1594         std::cout << "Nb nodes found : "<<nbFoundElems<<std::endl;
1595 #endif
1596         if (nbFoundElems ==0) {
1597           if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1598             newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
1599             anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1600           }
1601         }
1602         else if (nbFoundElems ==1) {
1603           const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1604           newId = (*aNodeToGhs3dIdMap.find(existingNode)).second;
1605           anExistingEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1606         }
1607         else
1608           isOK = false;
1609 #ifdef _DEBUG_
1610         std::cout << "GHS3D node ID: "<<newId<<std::endl;
1611 #endif
1612       }
1613       if (isOK)
1614         theKeptEnforcedEdges.insert(elem);
1615     }
1616   }
1617
1618   /* ENFORCED TRIANGLES ========================== */
1619
1620   // Iterate over the enforced triangles
1621   for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
1622     elem = elemIt->first;
1623     isOK = true;
1624     nodeIt = elem->nodesIterator();
1625     nbNodes = 3;
1626     while ( nodeIt->more() && nbNodes--) {
1627       // find GHS3D ID
1628       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1629       // Test if point is inside shape to mesh
1630       gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1631       TopAbs_State result = pntCls->GetPointState( myPoint );
1632       if ( result == TopAbs_OUT ) {
1633         isOK = false;
1634         break;
1635       }
1636       aNodeToTopAbs_StateMap.insert( make_pair( node, result ));
1637     }
1638     if (isOK) {
1639       nodeIt = elem->nodesIterator();
1640       nbNodes = 3;
1641       int newId = -1;
1642       while ( nodeIt->more() && nbNodes--) {
1643         // find GHS3D ID
1644         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1645         gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1646         nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1647 #ifdef _DEBUG_
1648         std::cout << "Nb nodes found : "<<nbFoundElems<<std::endl;
1649 #endif
1650         if (nbFoundElems ==0) {
1651           if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1652             newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
1653             anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1654           }
1655         }
1656         else if (nbFoundElems ==1) {
1657           const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1658           newId = (*aNodeToGhs3dIdMap.find(existingNode)).second;
1659           anExistingEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1660         }
1661         else
1662           isOK = false;
1663 #ifdef _DEBUG_
1664         std::cout << "GHS3D node ID: "<<newId<<std::endl;
1665 #endif
1666       }
1667       if (isOK)
1668         theKeptEnforcedTriangles.insert(elem);
1669     }
1670   }
1671
1672   // put nodes to theNodeByGhs3dId vector
1673 #ifdef _DEBUG_
1674   std::cout << "aNodeToGhs3dIdMap.size(): "<<aNodeToGhs3dIdMap.size()<<std::endl;
1675 #endif
1676   theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
1677   map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
1678   for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
1679   {
1680 //     std::cout << "n2id->first: "<<n2id->first<<std::endl;
1681     theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
1682   }
1683
1684   // put nodes to anEnforcedNodeToGhs3dIdMap vector
1685 #ifdef _DEBUG_
1686   std::cout << "anEnforcedNodeToGhs3dIdMap.size(): "<<anEnforcedNodeToGhs3dIdMap.size()<<std::endl;
1687 #endif
1688   theEnforcedNodeByGhs3dId.resize( anEnforcedNodeToGhs3dIdMap.size());
1689   n2id = anEnforcedNodeToGhs3dIdMap.begin();
1690   for ( ; n2id != anEnforcedNodeToGhs3dIdMap.end(); ++ n2id)
1691   {
1692     if (n2id->second > aNodeToGhs3dIdMap.size()) {
1693       theEnforcedNodeByGhs3dId[ n2id->second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // ghs3d ids count from 1
1694     }
1695   }
1696
1697
1698   /* ========================== NODES ========================== */
1699   vector<const SMDS_MeshNode*> theOrderedNodes, theRequiredNodes;
1700   std::set< std::vector<double> > nodesCoords;
1701   vector<const SMDS_MeshNode*>::const_iterator ghs3dNodeIt = theNodeByGhs3dId.begin();
1702   vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
1703
1704   (theNodeByGhs3dId.size() <= 1) ? tmpStr = " node" : " nodes";
1705   std::cout << theNodeByGhs3dId.size() << tmpStr << " from mesh ..." << std::endl;
1706   for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
1707   {
1708     const SMDS_MeshNode* node = *ghs3dNodeIt;
1709     std::vector<double> coords;
1710     coords.push_back(node->X());
1711     coords.push_back(node->Y());
1712     coords.push_back(node->Z());
1713     nodesCoords.insert(coords);
1714     theOrderedNodes.push_back(node);
1715   }
1716
1717   // Iterate over the enforced nodes given by enforced elements
1718   ghs3dNodeIt = theEnforcedNodeByGhs3dId.begin();
1719   after  = theEnforcedNodeByGhs3dId.end();
1720   (theEnforcedNodeByGhs3dId.size() <= 1) ? tmpStr = " node" : " nodes";
1721   std::cout << theEnforcedNodeByGhs3dId.size() << tmpStr << " from enforced elements ..." << std::endl;
1722   for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
1723   {
1724     const SMDS_MeshNode* node = *ghs3dNodeIt;
1725     std::vector<double> coords;
1726     coords.push_back(node->X());
1727     coords.push_back(node->Y());
1728     coords.push_back(node->Z());
1729 #ifdef _DEBUG_
1730     std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z();
1731 #endif
1732
1733     if (nodesCoords.find(coords) != nodesCoords.end()) {
1734       // node already exists in original mesh
1735 #ifdef _DEBUG_
1736       std::cout << " found" << std::endl;
1737 #endif
1738       continue;
1739     }
1740
1741     if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) {
1742       // node already exists in enforced vertices
1743 #ifdef _DEBUG_
1744       std::cout << " found" << std::endl;
1745 #endif
1746       continue;
1747     }
1748
1749 //     gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1750 //     nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1751 //     if (nbFoundElems ==0) {
1752 //       std::cout << " not found" << std::endl;
1753 //       if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1754 //         nodesCoords.insert(coords);
1755 //         theOrderedNodes.push_back(node);
1756 //       }
1757 //     }
1758 //     else {
1759 //       std::cout << " found in initial mesh" << std::endl;
1760 //       const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1761 //       nodesCoords.insert(coords);
1762 //       theOrderedNodes.push_back(existingNode);
1763 //     }
1764
1765 #ifdef _DEBUG_
1766     std::cout << " not found" << std::endl;
1767 #endif
1768
1769     nodesCoords.insert(coords);
1770     theOrderedNodes.push_back(node);
1771 //     theRequiredNodes.push_back(node);
1772   }
1773
1774
1775   // Iterate over the enforced nodes
1776   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator enfNodeIt;
1777   (theEnforcedNodes.size() <= 1) ? tmpStr = " node" : " nodes";
1778   std::cout << theEnforcedNodes.size() << tmpStr << " from enforced nodes ..." << std::endl;
1779   for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt)
1780   {
1781     const SMDS_MeshNode* node = enfNodeIt->first;
1782     std::vector<double> coords;
1783     coords.push_back(node->X());
1784     coords.push_back(node->Y());
1785     coords.push_back(node->Z());
1786 #ifdef _DEBUG_
1787     std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z();
1788 #endif
1789
1790     // Test if point is inside shape to mesh
1791     gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1792     TopAbs_State result = pntCls->GetPointState( myPoint );
1793     if ( result == TopAbs_OUT ) {
1794 #ifdef _DEBUG_
1795       std::cout << " out of volume" << std::endl;
1796 #endif
1797       continue;
1798     }
1799
1800     if (nodesCoords.find(coords) != nodesCoords.end()) {
1801 #ifdef _DEBUG_
1802       std::cout << " found in nodesCoords" << std::endl;
1803 #endif
1804 //       theRequiredNodes.push_back(node);
1805       continue;
1806     }
1807
1808     if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) {
1809 #ifdef _DEBUG_
1810       std::cout << " found in theEnforcedVertices" << std::endl;
1811 #endif
1812       continue;
1813     }
1814
1815 //     nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1816 //     if (nbFoundElems ==0) {
1817 //       std::cout << " not found" << std::endl;
1818 //       if (result == TopAbs_IN) {
1819 //         nodesCoords.insert(coords);
1820 //         theRequiredNodes.push_back(node);
1821 //       }
1822 //     }
1823 //     else {
1824 //       std::cout << " found in initial mesh" << std::endl;
1825 //       const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1826 // //       nodesCoords.insert(coords);
1827 //       theRequiredNodes.push_back(existingNode);
1828 //     }
1829 //
1830 //
1831 //
1832 //     if (pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems) == 0)
1833 //       continue;
1834
1835 //     if ( result != TopAbs_IN )
1836 //       continue;
1837
1838 #ifdef _DEBUG_
1839     std::cout << " not found" << std::endl;
1840 #endif
1841     nodesCoords.insert(coords);
1842 //     theOrderedNodes.push_back(node);
1843     theRequiredNodes.push_back(node);
1844   }
1845   int requiredNodes = theRequiredNodes.size();
1846
1847   int solSize = 0;
1848   std::vector<std::vector<double> > ReqVerTab;
1849   if (nbEnforcedVertices) {
1850 //    ReqVerTab.clear();
1851     (nbEnforcedVertices <= 1) ? tmpStr = " node" : " nodes";
1852     std::cout << nbEnforcedVertices << tmpStr << " from enforced vertices ..." << std::endl;
1853     // Iterate over the enforced vertices
1854     for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
1855       double x = vertexIt->first[0];
1856       double y = vertexIt->first[1];
1857       double z = vertexIt->first[2];
1858       // Test if point is inside shape to mesh
1859       gp_Pnt myPoint(x,y,z);
1860       TopAbs_State result = pntCls->GetPointState( myPoint );
1861       if ( result == TopAbs_OUT )
1862         continue;
1863       //if (pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems) == 0)
1864       //continue;
1865
1866 //       if ( result != TopAbs_IN )
1867 //         continue;
1868       std::vector<double> coords;
1869       coords.push_back(x);
1870       coords.push_back(y);
1871       coords.push_back(z);
1872       ReqVerTab.push_back(coords);
1873       enfVertexSizes.push_back(vertexIt->second);
1874       solSize++;
1875     }
1876   }
1877
1878
1879   // GmfVertices
1880   std::cout << "Begin writting required nodes in GmfVertices" << std::endl;
1881   std::cout << "Nb vertices: " << theOrderedNodes.size() << std::endl;
1882   GmfSetKwd(idx, GmfVertices, theOrderedNodes.size()/*+solSize*/);
1883   for (ghs3dNodeIt = theOrderedNodes.begin();ghs3dNodeIt != theOrderedNodes.end();++ghs3dNodeIt) {
1884     GmfSetLin(idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint);
1885   }
1886
1887   std::cout << "End writting required nodes in GmfVertices" << std::endl;
1888
1889   if (requiredNodes + solSize) {
1890     std::cout << "Begin writting in req and sol file" << std::endl;
1891     aNodeGroupByGhs3dId.resize( requiredNodes + solSize );
1892     idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1893     if (!idxRequired) {
1894       GmfCloseMesh(idx);
1895       return false;
1896     }
1897     idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1898     if (!idxSol) {
1899       GmfCloseMesh(idx);
1900       if (idxRequired)
1901         GmfCloseMesh(idxRequired);
1902       return false;
1903     }
1904     int TypTab[] = {GmfSca};
1905     double ValTab[] = {0.0};
1906     GmfSetKwd(idxRequired, GmfVertices, requiredNodes + solSize);
1907     GmfSetKwd(idxSol, GmfSolAtVertices, requiredNodes + solSize, 1, TypTab);
1908 //     int usedEnforcedNodes = 0;
1909 //     std::string gn = "";
1910     for (ghs3dNodeIt = theRequiredNodes.begin();ghs3dNodeIt != theRequiredNodes.end();++ghs3dNodeIt) {
1911       GmfSetLin(idxRequired, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint);
1912       GmfSetLin(idxSol, GmfSolAtVertices, ValTab);
1913       if (theEnforcedNodes.find((*ghs3dNodeIt)) != theEnforcedNodes.end())
1914         gn = theEnforcedNodes.find((*ghs3dNodeIt))->second;
1915       aNodeGroupByGhs3dId[usedEnforcedNodes] = gn;
1916       usedEnforcedNodes++;
1917     }
1918
1919     for (int i=0;i<solSize;i++) {
1920       std::cout << ReqVerTab[i][0] <<" "<< ReqVerTab[i][1] << " "<< ReqVerTab[i][2] << std::endl;
1921 #ifdef _DEBUG_
1922       std::cout << "enfVertexSizes.at("<<i<<"): " << enfVertexSizes.at(i) << std::endl;
1923 #endif
1924       double solTab[] = {enfVertexSizes.at(i)};
1925       GmfSetLin(idxRequired, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
1926       GmfSetLin(idxSol, GmfSolAtVertices, solTab);
1927       aNodeGroupByGhs3dId[usedEnforcedNodes] = enfVerticesWithGroup.find(ReqVerTab[i])->second;
1928 #ifdef _DEBUG_
1929       std::cout << "aNodeGroupByGhs3dId["<<usedEnforcedNodes<<"] = \""<<aNodeGroupByGhs3dId[usedEnforcedNodes]<<"\""<<std::endl;
1930 #endif
1931       usedEnforcedNodes++;
1932     }
1933     std::cout << "End writting in req and sol file" << std::endl;
1934   }
1935
1936   int nedge[2], ntri[3];
1937
1938   // GmfEdges
1939   int usedEnforcedEdges = 0;
1940   if (theKeptEnforcedEdges.size()) {
1941     anEdgeGroupByGhs3dId.resize( theKeptEnforcedEdges.size() );
1942 //    idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1943 //    if (!idxRequired)
1944 //      return false;
1945     GmfSetKwd(idx, GmfEdges, theKeptEnforcedEdges.size());
1946 //    GmfSetKwd(idxRequired, GmfEdges, theKeptEnforcedEdges.size());
1947     for(elemSetIt = theKeptEnforcedEdges.begin() ; elemSetIt != theKeptEnforcedEdges.end() ; ++elemSetIt) {
1948       elem = (*elemSetIt);
1949       nodeIt = elem->nodesIterator();
1950       int index=0;
1951       while ( nodeIt->more() ) {
1952         // find GHS3D ID
1953         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1954         map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
1955         if (it == anEnforcedNodeToGhs3dIdMap.end()) {
1956           it = anExistingEnforcedNodeToGhs3dIdMap.find(node);
1957           if (it == anEnforcedNodeToGhs3dIdMap.end())
1958             throw "Node not found";
1959         }
1960         nedge[index] = it->second;
1961         index++;
1962       }
1963       GmfSetLin(idx, GmfEdges, nedge[0], nedge[1], dummyint);
1964       anEdgeGroupByGhs3dId[usedEnforcedEdges] = theEnforcedEdges.find(elem)->second;
1965 //      GmfSetLin(idxRequired, GmfEdges, nedge[0], nedge[1], dummyint);
1966       usedEnforcedEdges++;
1967     }
1968 //    GmfCloseMesh(idxRequired);
1969   }
1970
1971
1972   if (usedEnforcedEdges) {
1973     GmfSetKwd(idx, GmfRequiredEdges, usedEnforcedEdges);
1974     for (int enfID=1;enfID<=usedEnforcedEdges;enfID++) {
1975       GmfSetLin(idx, GmfRequiredEdges, enfID);
1976     }
1977   }
1978
1979   // GmfTriangles
1980   int usedEnforcedTriangles = 0;
1981   if (anElemSet.size()+theKeptEnforcedTriangles.size()) {
1982     aFaceGroupByGhs3dId.resize( anElemSet.size()+theKeptEnforcedTriangles.size() );
1983     GmfSetKwd(idx, GmfTriangles, anElemSet.size()+theKeptEnforcedTriangles.size());
1984     int k=0;
1985     for(elemSetIt = anElemSet.begin() ; elemSetIt != anElemSet.end() ; ++elemSetIt,++k) {
1986       elem = (*elemSetIt);
1987       nodeIt = elem->nodesIterator();
1988       int index=0;
1989       for ( int j = 0; j < 3; ++j ) {
1990         // find GHS3D ID
1991         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1992         map< const SMDS_MeshNode*,int >::iterator it = aNodeToGhs3dIdMap.find(node);
1993         if (it == aNodeToGhs3dIdMap.end())
1994           throw "Node not found";
1995         ntri[index] = it->second;
1996         index++;
1997       }
1998       GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint);
1999       aFaceGroupByGhs3dId[k] = "";
2000     }
2001     if (theKeptEnforcedTriangles.size()) {
2002       for(elemSetIt = theKeptEnforcedTriangles.begin() ; elemSetIt != theKeptEnforcedTriangles.end() ; ++elemSetIt,++k) {
2003         elem = (*elemSetIt);
2004         nodeIt = elem->nodesIterator();
2005         int index=0;
2006         for ( int j = 0; j < 3; ++j ) {
2007           // find GHS3D ID
2008           const SMDS_MeshNode* node = castToNode( nodeIt->next() );
2009           map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
2010           if (it == anEnforcedNodeToGhs3dIdMap.end()) {
2011             it = anExistingEnforcedNodeToGhs3dIdMap.find(node);
2012             if (it == anEnforcedNodeToGhs3dIdMap.end())
2013               throw "Node not found";
2014           }
2015           ntri[index] = it->second;
2016           index++;
2017         }
2018         GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint);
2019         aFaceGroupByGhs3dId[k] = theEnforcedTriangles.find(elem)->second;
2020         usedEnforcedTriangles++;
2021       }
2022     }
2023   }
2024
2025
2026   if (usedEnforcedTriangles) {
2027     GmfSetKwd(idx, GmfRequiredTriangles, usedEnforcedTriangles);
2028     for (int enfID=1;enfID<=usedEnforcedTriangles;enfID++)
2029       GmfSetLin(idx, GmfRequiredTriangles, anElemSet.size()+enfID);
2030   }
2031
2032   GmfCloseMesh(idx);
2033   if (idxRequired)
2034     GmfCloseMesh(idxRequired);
2035   if (idxSol)
2036     GmfCloseMesh(idxSol);
2037
2038   return true;
2039
2040 }
2041
2042 // static bool writeGMFFile(const char*                                    theMeshFileName,
2043 //                         const char*                                     theRequiredFileName,
2044 //                         const char*                                     theSolFileName,
2045 //                         SMESH_MesherHelper&                             theHelper,
2046 //                         const SMESH_ProxyMesh&                          theProxyMesh,
2047 //                         std::map <int,int> &                            theNodeId2NodeIndexMap,
2048 //                         std::map <int,int> &                            theSmdsToGhs3dIdMap,
2049 //                         std::map <int,const SMDS_MeshNode*> &           theGhs3dIdToNodeMap,
2050 //                         TIDSortedNodeSet &                              theEnforcedNodes,
2051 //                         TIDSortedElemSet &                              theEnforcedEdges,
2052 //                         TIDSortedElemSet &                              theEnforcedTriangles,
2053 // //                         TIDSortedElemSet &                              theEnforcedQuadrangles,
2054 //                         GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices)
2055 // {
2056 //   MESSAGE("writeGMFFile with geometry");
2057 //   int idx, idxRequired, idxSol;
2058 //   int nbv, nbev, nben, aGhs3dID = 0;
2059 //   const int dummyint = 0;
2060 //   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt;
2061 //   std::vector<double> enfVertexSizes;
2062 //   TIDSortedNodeSet::const_iterator enfNodeIt;
2063 //   const SMDS_MeshNode* node;
2064 //   SMDS_NodeIteratorPtr nodeIt;
2065 //
2066 //   idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
2067 //   if (!idx)
2068 //     return false;
2069 //
2070 //   SMESHDS_Mesh * theMeshDS = theHelper.GetMeshDS();
2071 //
2072 //   /* ========================== NODES ========================== */
2073 //   // NB_NODES
2074 //   nbv = theMeshDS->NbNodes();
2075 //   if ( nbv == 0 )
2076 //     return false;
2077 //   nbev = theEnforcedVertices.size();
2078 //   nben = theEnforcedNodes.size();
2079 //
2080 //   // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
2081 //   // The problem is in nodes on degenerated edges, we need to skip nodes which are free
2082 //   // and replace not-free nodes on edges by the node on vertex
2083 //   TNodeNodeMap n2nDegen; // map a node on degenerated edge to a node on vertex
2084 //   TNodeNodeMap::iterator n2nDegenIt;
2085 //   if ( theHelper.HasDegeneratedEdges() )
2086 //   {
2087 //     set<int> checkedSM;
2088 //     for (TopExp_Explorer e(theMeshDS->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
2089 //     {
2090 //       SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
2091 //       if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() ))
2092 //       {
2093 //         if ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() )
2094 //         {
2095 //           TopoDS_Shape vertex = TopoDS_Iterator( e.Current() ).Value();
2096 //           const SMDS_MeshNode* vNode = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), theMeshDS);
2097 //           {
2098 //             SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2099 //             while ( nIt->more() )
2100 //               n2nDegen.insert( make_pair( nIt->next(), vNode ));
2101 //           }
2102 //         }
2103 //       }
2104 //     }
2105 //   }
2106 //
2107 //   const bool isQuadMesh =
2108 //     theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
2109 //     theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
2110 //     theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC );
2111 //
2112 //   std::vector<std::vector<double> > VerTab;
2113 //   std::set<std::vector<double> > VerMap;
2114 //   VerTab.clear();
2115 //   std::vector<double> aVerTab;
2116 //   // Loop from 1 to NB_NODES
2117 //
2118 //   nodeIt = theMeshDS->nodesIterator();
2119 //
2120 //   while ( nodeIt->more() )
2121 //   {
2122 //     node = nodeIt->next();
2123 //     if ( isQuadMesh && theHelper.IsMedium( node )) // Issue 0021238
2124 //       continue;
2125 //     if ( n2nDegen.count( node ) ) // Issue 0020674
2126 //       continue;
2127 //
2128 //     std::vector<double> coords;
2129 //     coords.push_back(node->X());
2130 //     coords.push_back(node->Y());
2131 //     coords.push_back(node->Z());
2132 //     if (VerMap.find(coords) != VerMap.end()) {
2133 //       aGhs3dID = theSmdsToGhs3dIdMap[node->GetID()];
2134 //       theGhs3dIdToNodeMap[theSmdsToGhs3dIdMap[node->GetID()]] = node;
2135 //       continue;
2136 //     }
2137 //     VerTab.push_back(coords);
2138 //     VerMap.insert(coords);
2139 //     aGhs3dID++;
2140 //     theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
2141 //     theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
2142 //   }
2143 //
2144 //
2145 //   /* ENFORCED NODES ========================== */
2146 //   if (nben) {
2147 //     std::cout << "Add " << nben << " enforced nodes to input .mesh file" << std::endl;
2148 //     for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt) {
2149 //       double x = (*enfNodeIt)->X();
2150 //       double y = (*enfNodeIt)->Y();
2151 //       double z = (*enfNodeIt)->Z();
2152 //       // Test if point is inside shape to mesh
2153 //       gp_Pnt myPoint(x,y,z);
2154 //       BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh());
2155 //       scl.Perform(myPoint, 1e-7);
2156 //       TopAbs_State result = scl.State();
2157 //       if ( result != TopAbs_IN )
2158 //         continue;
2159 //       std::vector<double> coords;
2160 //       coords.push_back(x);
2161 //       coords.push_back(y);
2162 //       coords.push_back(z);
2163 //       if (theEnforcedVertices.find(coords) != theEnforcedVertices.end())
2164 //         continue;
2165 //       if (VerMap.find(coords) != VerMap.end())
2166 //         continue;
2167 //       VerTab.push_back(coords);
2168 //       VerMap.insert(coords);
2169 //       aGhs3dID++;
2170 //       theNodeId2NodeIndexMap.insert( make_pair( (*enfNodeIt)->GetID(), aGhs3dID ));
2171 //     }
2172 //   }
2173 //
2174 //
2175 //   /* ENFORCED VERTICES ========================== */
2176 //   int solSize = 0;
2177 //   std::vector<std::vector<double> > ReqVerTab;
2178 //   ReqVerTab.clear();
2179 //   if (nbev) {
2180 //     std::cout << "Add " << nbev << " enforced vertices to input .mesh file" << std::endl;
2181 //     for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
2182 //       double x = vertexIt->first[0];
2183 //       double y = vertexIt->first[1];
2184 //       double z = vertexIt->first[2];
2185 //       // Test if point is inside shape to mesh
2186 //       gp_Pnt myPoint(x,y,z);
2187 //       BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh());
2188 //       scl.Perform(myPoint, 1e-7);
2189 //       TopAbs_State result = scl.State();
2190 //       if ( result != TopAbs_IN )
2191 //         continue;
2192 //       enfVertexSizes.push_back(vertexIt->second);
2193 //       std::vector<double> coords;
2194 //       coords.push_back(x);
2195 //       coords.push_back(y);
2196 //       coords.push_back(z);
2197 //       if (VerMap.find(coords) != VerMap.end())
2198 //         continue;
2199 //       ReqVerTab.push_back(coords);
2200 //       VerMap.insert(coords);
2201 //       solSize++;
2202 //     }
2203 //   }
2204 //
2205 //
2206 //   /* ========================== FACES ========================== */
2207 //
2208 //   int nbTriangles = 0/*, nbQuadrangles = 0*/, aSmdsID;
2209 //   TopTools_IndexedMapOfShape facesMap, trianglesMap/*, quadranglesMap*/;
2210 //   TIDSortedElemSet::const_iterator elemIt;
2211 //   const SMESHDS_SubMesh* theSubMesh;
2212 //   TopoDS_Shape aShape;
2213 //   SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
2214 //   const SMDS_MeshElement* aFace;
2215 //   map<int,int>::const_iterator itOnMap;
2216 //   std::vector<std::vector<int> > tt, qt,et;
2217 //   tt.clear();
2218 //   qt.clear();
2219 //   et.clear();
2220 //   std::vector<int> att, aqt, aet;
2221 //
2222 //   TopExp::MapShapes( theMeshDS->ShapeToMesh(), TopAbs_FACE, facesMap );
2223 //
2224 //   for ( int i = 1; i <= facesMap.Extent(); ++i )
2225 //     if (( theSubMesh  = theProxyMesh.GetSubMesh( facesMap(i))))
2226 //     {
2227 //       SMDS_ElemIteratorPtr it = theSubMesh->GetElements();
2228 //       while (it->more())
2229 //       {
2230 //         const SMDS_MeshElement *elem = it->next();
2231 //         int nbCornerNodes = elem->NbCornerNodes();
2232 //         if (nbCornerNodes == 3)
2233 //         {
2234 //           trianglesMap.Add(facesMap(i));
2235 //           nbTriangles ++;
2236 //         }
2237 // //         else if (nbCornerNodes == 4)
2238 // //         {
2239 // //           quadranglesMap.Add(facesMap(i));
2240 // //           nbQuadrangles ++;
2241 // //         }
2242 //       }
2243 //     }
2244 //
2245 //   /* TRIANGLES ========================== */
2246 //   if (nbTriangles) {
2247 //     for ( int i = 1; i <= trianglesMap.Extent(); i++ )
2248 //     {
2249 //       aShape = trianglesMap(i);
2250 //       theSubMesh = theProxyMesh.GetSubMesh(aShape);
2251 //       if ( !theSubMesh ) continue;
2252 //       itOnSubMesh = theSubMesh->GetElements();
2253 //       while ( itOnSubMesh->more() )
2254 //       {
2255 //         aFace = itOnSubMesh->next();
2256 //         itOnSubFace = aFace->nodesIterator();
2257 //         att.clear();
2258 //         for ( int j = 0; j < 3; ++j ) {
2259 //           // find GHS3D ID
2260 //           node = castToNode( itOnSubFace->next() );
2261 //           if (( n2nDegenIt = n2nDegen.find( node )) != n2nDegen.end() )
2262 //             node = n2nDegenIt->second;
2263 //           aSmdsID = node->GetID();
2264 //           itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
2265 //           ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
2266 //           att.push_back((*itOnMap).second);
2267 //         }
2268 //         tt.push_back(att);
2269 //       }
2270 //     }
2271 //   }
2272 //
2273 //   if (theEnforcedTriangles.size()) {
2274 //     std::cout << "Add " << theEnforcedTriangles.size() << " enforced triangles to input .mesh file" << std::endl;
2275 //     // Iterate over the enforced triangles
2276 //     for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
2277 //       aFace = (*elemIt);
2278 //       itOnSubFace = aFace->nodesIterator();
2279 //       bool isOK = true;
2280 //       att.clear();
2281 //
2282 //       for ( int j = 0; j < 3; ++j ) {
2283 //         node = castToNode( itOnSubFace->next() );
2284 //         if (( n2nDegenIt = n2nDegen.find( node )) != n2nDegen.end() )
2285 //           node = n2nDegenIt->second;
2286 // //         std::cout << node;
2287 //         double x = node->X();
2288 //         double y = node->Y();
2289 //         double z = node->Z();
2290 //         // Test if point is inside shape to mesh
2291 //         gp_Pnt myPoint(x,y,z);
2292 //         BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh());
2293 //         scl.Perform(myPoint, 1e-7);
2294 //         TopAbs_State result = scl.State();
2295 //         if ( result != TopAbs_IN ) {
2296 //           isOK = false;
2297 //           theEnforcedTriangles.erase(elemIt);
2298 //           continue;
2299 //         }
2300 //         std::vector<double> coords;
2301 //         coords.push_back(x);
2302 //         coords.push_back(y);
2303 //         coords.push_back(z);
2304 //         if (VerMap.find(coords) != VerMap.end()) {
2305 //           att.push_back(theNodeId2NodeIndexMap[node->GetID()]);
2306 //           continue;
2307 //         }
2308 //         VerTab.push_back(coords);
2309 //         VerMap.insert(coords);
2310 //         aGhs3dID++;
2311 //         theNodeId2NodeIndexMap.insert( make_pair( node->GetID(), aGhs3dID ));
2312 //         att.push_back(aGhs3dID);
2313 //       }
2314 //       if (isOK)
2315 //         tt.push_back(att);
2316 //     }
2317 //   }
2318 //
2319 //
2320 //   /* ========================== EDGES ========================== */
2321 //
2322 //   if (theEnforcedEdges.size()) {
2323 //     // Iterate over the enforced edges
2324 //     std::cout << "Add " << theEnforcedEdges.size() << " enforced edges to input .mesh file" << std::endl;
2325 //     for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
2326 //       aFace = (*elemIt);
2327 //       bool isOK = true;
2328 //       itOnSubFace = aFace->nodesIterator();
2329 //       aet.clear();
2330 //       for ( int j = 0; j < 2; ++j ) {
2331 //         node = castToNode( itOnSubFace->next() );
2332 //         if (( n2nDegenIt = n2nDegen.find( node )) != n2nDegen.end() )
2333 //           node = n2nDegenIt->second;
2334 //         double x = node->X();
2335 //         double y = node->Y();
2336 //         double z = node->Z();
2337 //         // Test if point is inside shape to mesh
2338 //         gp_Pnt myPoint(x,y,z);
2339 //         BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh());
2340 //         scl.Perform(myPoint, 1e-7);
2341 //         TopAbs_State result = scl.State();
2342 //         if ( result != TopAbs_IN ) {
2343 //           isOK = false;
2344 //           theEnforcedEdges.erase(elemIt);
2345 //           continue;
2346 //         }
2347 //         std::vector<double> coords;
2348 //         coords.push_back(x);
2349 //         coords.push_back(y);
2350 //         coords.push_back(z);
2351 //         if (VerMap.find(coords) != VerMap.end()) {
2352 //           aet.push_back(theNodeId2NodeIndexMap[node->GetID()]);
2353 //           continue;
2354 //         }
2355 //         VerTab.push_back(coords);
2356 //         VerMap.insert(coords);
2357 //
2358 //         aGhs3dID++;
2359 //         theNodeId2NodeIndexMap.insert( make_pair( node->GetID(), aGhs3dID ));
2360 //         aet.push_back(aGhs3dID);
2361 //       }
2362 //       if (isOK)
2363 //         et.push_back(aet);
2364 //     }
2365 //   }
2366 //
2367 //
2368 //   /* Write vertices number */
2369 //   MESSAGE("Number of vertices: "<<aGhs3dID);
2370 //   MESSAGE("Size of vector: "<<VerTab.size());
2371 //   GmfSetKwd(idx, GmfVertices, aGhs3dID/*+solSize*/);
2372 //   for (int i=0;i<aGhs3dID;i++)
2373 //     GmfSetLin(idx, GmfVertices, VerTab[i][0], VerTab[i][1], VerTab[i][2], dummyint);
2374 // //   for (int i=0;i<solSize;i++) {
2375 // //     std::cout << ReqVerTab[i][0] <<" "<< ReqVerTab[i][1] << " "<< ReqVerTab[i][2] << std::endl;
2376 // //     GmfSetLin(idx, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
2377 // //   }
2378 //
2379 //   if (solSize) {
2380 //     idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
2381 //     if (!idxRequired) {
2382 //       GmfCloseMesh(idx);
2383 //       return false;
2384 //     }
2385 //     idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
2386 //     if (!idxSol){
2387 //       GmfCloseMesh(idx);
2388 //       if (idxRequired)
2389 //         GmfCloseMesh(idxRequired);
2390 //       return false;
2391 //     }
2392 //
2393 //     int TypTab[] = {GmfSca};
2394 //     GmfSetKwd(idxRequired, GmfVertices, solSize);
2395 //     GmfSetKwd(idxSol, GmfSolAtVertices, solSize, 1, TypTab);
2396 //
2397 //     for (int i=0;i<solSize;i++) {
2398 //       double solTab[] = {enfVertexSizes.at(i)};
2399 //       GmfSetLin(idxRequired, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
2400 //       GmfSetLin(idxSol, GmfSolAtVertices, solTab);
2401 //     }
2402 //     GmfCloseMesh(idxRequired);
2403 //     GmfCloseMesh(idxSol);
2404 //   }
2405 //
2406 //   /* Write triangles number */
2407 //   if (tt.size()) {
2408 //     GmfSetKwd(idx, GmfTriangles, tt.size());
2409 //     for (int i=0;i<tt.size();i++)
2410 //       GmfSetLin(idx, GmfTriangles, tt[i][0], tt[i][1], tt[i][2], dummyint);
2411 //   }
2412 //
2413 //   /* Write edges number */
2414 //   if (et.size()) {
2415 //     GmfSetKwd(idx, GmfEdges, et.size());
2416 //     for (int i=0;i<et.size();i++)
2417 //       GmfSetLin(idx, GmfEdges, et[i][0], et[i][1], dummyint);
2418 //   }
2419 //
2420 //   /* QUADRANGLES ========================== */
2421 //   // TODO: add pyramids ?
2422 // //   if (nbQuadrangles) {
2423 // //     for ( int i = 1; i <= quadranglesMap.Extent(); i++ )
2424 // //     {
2425 // //       aShape = quadranglesMap(i);
2426 // //       theSubMesh = theProxyMesh.GetSubMesh(aShape);
2427 // //       if ( !theSubMesh ) continue;
2428 // //       itOnSubMesh = theSubMesh->GetElements();
2429 // //       for ( int j = 0; j < 4; ++j )
2430 // //       {
2431 // //         aFace = itOnSubMesh->next();
2432 // //         itOnSubFace = aFace->nodesIterator();
2433 // //         aqt.clear();
2434 // //         while ( itOnSubFace->more() ) {
2435 // //           // find GHS3D ID
2436 // //           aSmdsID = itOnSubFace->next()->GetID();
2437 // //           itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
2438 // //           ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
2439 // //           aqt.push_back((*itOnMap).second);
2440 // //         }
2441 // //         qt.push_back(aqt);
2442 // //       }
2443 // //     }
2444 // //   }
2445 // //
2446 // //   if (theEnforcedQuadrangles.size()) {
2447 // //     // Iterate over the enforced triangles
2448 // //     for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) {
2449 // //       aFace = (*elemIt);
2450 // //       bool isOK = true;
2451 // //       itOnSubFace = aFace->nodesIterator();
2452 // //       aqt.clear();
2453 // //       for ( int j = 0; j < 4; ++j ) {
2454 // //         int aNodeID = itOnSubFace->next()->GetID();
2455 // //         itOnMap = theNodeId2NodeIndexMap.find(aNodeID);
2456 // //         if (itOnMap != theNodeId2NodeIndexMap.end())
2457 // //           aqt.push_back((*itOnMap).second);
2458 // //         else {
2459 // //           isOK = false;
2460 // //           theEnforcedQuadrangles.erase(elemIt);
2461 // //           break;
2462 // //         }
2463 // //       }
2464 // //       if (isOK)
2465 // //         qt.push_back(aqt);
2466 // //     }
2467 // //   }
2468 // //
2469 //
2470 // //   /* Write quadrilaterals number */
2471 // //   if (qt.size()) {
2472 // //     GmfSetKwd(idx, GmfQuadrilaterals, qt.size());
2473 // //     for (int i=0;i<qt.size();i++)
2474 // //       GmfSetLin(idx, GmfQuadrilaterals, qt[i][0], qt[i][1], qt[i][2], qt[i][3], dummyint);
2475 // //   }
2476 //
2477 //   GmfCloseMesh(idx);
2478 //   return true;
2479 // }
2480
2481
2482 //=======================================================================
2483 //function : writeFaces
2484 //purpose  :
2485 //=======================================================================
2486
2487 static bool writeFaces (ofstream &              theFile,
2488                         const SMESH_ProxyMesh&  theMesh,
2489                         const TopoDS_Shape&     theShape,
2490                         const map <int,int> &   theSmdsToGhs3dIdMap,
2491                         const map <int,int> &   theEnforcedNodeIdToGhs3dIdMap,
2492                         GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges,
2493                         GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles)
2494 {
2495   // record structure:
2496   //
2497   // NB_ELEMS DUMMY_INT
2498   // Loop from 1 to NB_ELEMS
2499   // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
2500
2501   TopoDS_Shape aShape;
2502   const SMESHDS_SubMesh* theSubMesh;
2503   const SMDS_MeshElement* aFace;
2504   const char* space    = "  ";
2505   const int   dummyint = 0;
2506   map<int,int>::const_iterator itOnMap;
2507   SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
2508   int nbNodes, aSmdsID;
2509
2510   TIDSortedElemSet::const_iterator elemIt;
2511   int nbEnforcedEdges       = theEnforcedEdges.size();
2512   int nbEnforcedTriangles   = theEnforcedTriangles.size();
2513
2514   // count triangles bound to geometry
2515   int nbTriangles = 0;
2516
2517   TopTools_IndexedMapOfShape facesMap, trianglesMap;
2518   TopExp::MapShapes( theShape, TopAbs_FACE, facesMap );
2519
2520   int nbFaces = facesMap.Extent();
2521
2522   for ( int i = 1; i <= nbFaces; ++i )
2523     if (( theSubMesh  = theMesh.GetSubMesh( facesMap(i))))
2524       nbTriangles += theSubMesh->NbElements();
2525   std::string tmpStr;
2526   (nbFaces == 0 || nbFaces == 1) ? tmpStr = " shape " : tmpStr = " shapes " ;
2527   std::cout << "    " << nbFaces << tmpStr << "of 2D dimension";
2528   int nbEnforcedElements = nbEnforcedEdges+nbEnforcedTriangles;
2529   if (nbEnforcedElements > 0) {
2530     (nbEnforcedElements == 1) ? tmpStr = "shape:" : tmpStr = "shapes:";
2531     std::cout << " and" << std::endl;
2532     std::cout << "    " << nbEnforcedElements
2533                         << " enforced " << tmpStr << std::endl;
2534   }
2535   else
2536     std::cout << std::endl;
2537   if (nbEnforcedEdges) {
2538     (nbEnforcedEdges == 1) ? tmpStr = "edge" : tmpStr = "edges";
2539     std::cout << "      " << nbEnforcedEdges << " enforced " << tmpStr << std::endl;
2540   }
2541   if (nbEnforcedTriangles) {
2542     (nbEnforcedTriangles == 1) ? tmpStr = "triangle" : tmpStr = "triangles";
2543     std::cout << "      " << nbEnforcedTriangles << " enforced " << tmpStr << std::endl;
2544   }
2545   std::cout << std::endl;
2546
2547 //   theFile << space << nbTriangles << space << dummyint << std::endl;
2548   std::ostringstream globalStream, localStream, aStream;
2549
2550   for ( int i = 1; i <= facesMap.Extent(); i++ )
2551   {
2552     aShape = facesMap(i);
2553     theSubMesh = theMesh.GetSubMesh(aShape);
2554     if ( !theSubMesh ) continue;
2555     itOnSubMesh = theSubMesh->GetElements();
2556     while ( itOnSubMesh->more() )
2557     {
2558       aFace = itOnSubMesh->next();
2559       nbNodes = aFace->NbCornerNodes();
2560
2561       localStream << nbNodes << space;
2562
2563       itOnSubFace = aFace->nodesIterator();
2564       for ( int j = 0; j < 3; ++j ) {
2565         // find GHS3D ID
2566         aSmdsID = itOnSubFace->next()->GetID();
2567         itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
2568         // if ( itOnMap == theSmdsToGhs3dIdMap.end() ) {
2569         //   cout << "not found node: " << aSmdsID << endl;
2570         //   return false;
2571         // }
2572         ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
2573
2574         localStream << (*itOnMap).second << space ;
2575       }
2576
2577       // (NB_NODES + 1) times: DUMMY_INT
2578       for ( int j=0; j<=nbNodes; j++)
2579         localStream << dummyint << space ;
2580
2581       localStream << std::endl;
2582     }
2583   }
2584
2585   globalStream << localStream.str();
2586   localStream.str("");
2587
2588   //
2589   //        FACES : END
2590   //
2591
2592 //   //
2593 //   //        ENFORCED EDGES : BEGIN
2594 //   //
2595 //
2596 //   // Iterate over the enforced edges
2597 //   int usedEnforcedEdges = 0;
2598 //   bool isOK;
2599 //   for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
2600 //     aFace = (*elemIt);
2601 //     isOK = true;
2602 //     itOnSubFace = aFace->nodesIterator();
2603 //     aStream.str("");
2604 //     aStream << "2" << space ;
2605 //     for ( int j = 0; j < 2; ++j ) {
2606 //       aSmdsID = itOnSubFace->next()->GetID();
2607 //       itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID);
2608 //       if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end())
2609 //         aStream << (*itOnMap).second << space;
2610 //       else {
2611 //         isOK = false;
2612 //         break;
2613 //       }
2614 //     }
2615 //     if (isOK) {
2616 //       for ( int j=0; j<=2; j++)
2617 //         aStream << dummyint << space ;
2618 // //       aStream << dummyint << space << dummyint;
2619 //       localStream << aStream.str() << std::endl;
2620 //       usedEnforcedEdges++;
2621 //     }
2622 //   }
2623 //
2624 //   if (usedEnforcedEdges) {
2625 //     globalStream << localStream.str();
2626 //     localStream.str("");
2627 //   }
2628 //
2629 //   //
2630 //   //        ENFORCED EDGES : END
2631 //   //
2632 //   //
2633 //
2634 //   //
2635 //   //        ENFORCED TRIANGLES : BEGIN
2636 //   //
2637 //     // Iterate over the enforced triangles
2638 //   int usedEnforcedTriangles = 0;
2639 //   for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
2640 //     aFace = (*elemIt);
2641 //     nbNodes = aFace->NbCornerNodes();
2642 //     isOK = true;
2643 //     itOnSubFace = aFace->nodesIterator();
2644 //     aStream.str("");
2645 //     aStream << nbNodes << space ;
2646 //     for ( int j = 0; j < 3; ++j ) {
2647 //       aSmdsID = itOnSubFace->next()->GetID();
2648 //       itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID);
2649 //       if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end())
2650 //         aStream << (*itOnMap).second << space;
2651 //       else {
2652 //         isOK = false;
2653 //         break;
2654 //       }
2655 //     }
2656 //     if (isOK) {
2657 //       for ( int j=0; j<=3; j++)
2658 //         aStream << dummyint << space ;
2659 //       localStream << aStream.str() << std::endl;
2660 //       usedEnforcedTriangles++;
2661 //     }
2662 //   }
2663 //
2664 //   if (usedEnforcedTriangles) {
2665 //     globalStream << localStream.str();
2666 //     localStream.str("");
2667 //   }
2668 //
2669 //   //
2670 //   //        ENFORCED TRIANGLES : END
2671 //   //
2672
2673   theFile
2674   << nbTriangles/*+usedEnforcedTriangles+usedEnforcedEdges*/
2675   << " 0" << std::endl
2676   << globalStream.str();
2677
2678   return true;
2679 }
2680
2681 //=======================================================================
2682 //function : writePoints
2683 //purpose  :
2684 //=======================================================================
2685
2686 static bool writePoints (ofstream &                       theFile,
2687                          SMESH_MesherHelper&              theHelper,
2688                          map <int,int> &                  theSmdsToGhs3dIdMap,
2689                          map <int,int> &                  theEnforcedNodeIdToGhs3dIdMap,
2690                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap,
2691                          GHS3DPlugin_Hypothesis::TID2SizeMap & theNodeIDToSizeMap,
2692                          GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices,
2693                          GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap & theEnforcedNodes,
2694                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges,
2695                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles)
2696 {
2697   // record structure:
2698   //
2699   // NB_NODES
2700   // Loop from 1 to NB_NODES
2701   //   X Y Z DUMMY_INT
2702
2703   SMESHDS_Mesh * theMeshDS = theHelper.GetMeshDS();
2704   int nbNodes = theMeshDS->NbNodes();
2705   if ( nbNodes == 0 )
2706     return false;
2707
2708   int nbEnforcedVertices = theEnforcedVertices.size();
2709   int nbEnforcedNodes    = theEnforcedNodes.size();
2710
2711   const TopoDS_Shape shapeToMesh = theMeshDS->ShapeToMesh();
2712
2713   int aGhs3dID = 1;
2714   SMDS_NodeIteratorPtr nodeIt = theMeshDS->nodesIterator();
2715   const SMDS_MeshNode* node;
2716
2717   // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
2718   // The problem is in nodes on degenerated edges, we need to skip nodes which are free
2719   // and replace not-free nodes on degenerated edges by the node on vertex
2720   TNodeNodeMap n2nDegen; // map a node on degenerated edge to a node on vertex
2721   TNodeNodeMap::iterator n2nDegenIt;
2722   if ( theHelper.HasDegeneratedEdges() )
2723   {
2724     set<int> checkedSM;
2725     for (TopExp_Explorer e(theMeshDS->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
2726     {
2727       SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
2728       if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() ))
2729       {
2730         if ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() )
2731         {
2732           TopoDS_Shape vertex = TopoDS_Iterator( e.Current() ).Value();
2733           const SMDS_MeshNode* vNode = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), theMeshDS);
2734           {
2735             SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2736             while ( nIt->more() )
2737               n2nDegen.insert( make_pair( nIt->next(), vNode ));
2738           }
2739         }
2740       }
2741     }
2742     nbNodes -= n2nDegen.size();
2743   }
2744
2745   const bool isQuadMesh =
2746     theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
2747     theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
2748     theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC );
2749   if ( isQuadMesh )
2750   {
2751     // descrease nbNodes by nb of medium nodes
2752     while ( nodeIt->more() )
2753     {
2754       node = nodeIt->next();
2755       if ( !theHelper.IsDegenShape( node->getshapeId() ))
2756         nbNodes -= int( theHelper.IsMedium( node ));
2757     }
2758     nodeIt = theMeshDS->nodesIterator();
2759   }
2760
2761   const char* space    = "  ";
2762   const int   dummyint = 0;
2763
2764   std::string tmpStr;
2765   (nbNodes == 0 || nbNodes == 1) ? tmpStr = " node" : tmpStr = " nodes";
2766   // NB_NODES
2767   std::cout << std::endl;
2768   std::cout << "The initial 2D mesh contains :" << std::endl;
2769   std::cout << "    " << nbNodes << tmpStr << std::endl;
2770   if (nbEnforcedVertices > 0) {
2771     (nbEnforcedVertices == 1) ? tmpStr = "vertex" : tmpStr = "vertices";
2772     std::cout << "    " << nbEnforcedVertices << " enforced " << tmpStr << std::endl;
2773   }
2774   if (nbEnforcedNodes > 0) {
2775     (nbEnforcedNodes == 1) ? tmpStr = "node" : tmpStr = "nodes";
2776     std::cout << "    " << nbEnforcedNodes << " enforced " << tmpStr << std::endl;
2777   }
2778   std::cout << std::endl;
2779   std::cout << "Start writing in 'points' file ..." << std::endl;
2780
2781   theFile << nbNodes << std::endl;
2782
2783   // Loop from 1 to NB_NODES
2784
2785   while ( nodeIt->more() )
2786   {
2787     node = nodeIt->next();
2788     if ( isQuadMesh && theHelper.IsMedium( node )) // Issue 0021238
2789       continue;
2790     if ( n2nDegen.count( node ) ) // Issue 0020674
2791       continue;
2792
2793     theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
2794     theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
2795     aGhs3dID++;
2796
2797     // X Y Z DUMMY_INT
2798     theFile
2799     << node->X() << space
2800     << node->Y() << space
2801     << node->Z() << space
2802     << dummyint;
2803
2804     theFile << std::endl;
2805
2806   }
2807
2808   // Iterate over the enforced nodes
2809   std::map<int,double> enfVertexIndexSizeMap;
2810   if (nbEnforcedNodes) {
2811     GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator nodeIt = theEnforcedNodes.begin();
2812     for( ; nodeIt != theEnforcedNodes.end() ; ++nodeIt) {
2813       double x = nodeIt->first->X();
2814       double y = nodeIt->first->Y();
2815       double z = nodeIt->first->Z();
2816       // Test if point is inside shape to mesh
2817       gp_Pnt myPoint(x,y,z);
2818       BRepClass3d_SolidClassifier scl(shapeToMesh);
2819       scl.Perform(myPoint, 1e-7);
2820       TopAbs_State result = scl.State();
2821       if ( result != TopAbs_IN )
2822         continue;
2823       std::vector<double> coords;
2824       coords.push_back(x);
2825       coords.push_back(y);
2826       coords.push_back(z);
2827       if (theEnforcedVertices.find(coords) != theEnforcedVertices.end())
2828         continue;
2829
2830 //      double size = theNodeIDToSizeMap.find(nodeIt->first->GetID())->second;
2831   //       theGhs3dIdToNodeMap.insert( make_pair( nbNodes + i, (*nodeIt) ));
2832   //       MESSAGE("Adding enforced node (" << x << "," << y <<"," << z << ")");
2833       // X Y Z PHY_SIZE DUMMY_INT
2834       theFile
2835       << x << space
2836       << y << space
2837       << z << space
2838       << -1 << space
2839       << dummyint << space;
2840       theFile << std::endl;
2841       theEnforcedNodeIdToGhs3dIdMap.insert( make_pair( nodeIt->first->GetID(), aGhs3dID ));
2842       enfVertexIndexSizeMap[aGhs3dID] = -1;
2843       aGhs3dID++;
2844   //     else
2845   //         MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added ");
2846     }
2847   }
2848
2849   if (nbEnforcedVertices) {
2850     // Iterate over the enforced vertices
2851     GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt = theEnforcedVertices.begin();
2852     for( ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
2853       double x = vertexIt->first[0];
2854       double y = vertexIt->first[1];
2855       double z = vertexIt->first[2];
2856       // Test if point is inside shape to mesh
2857       gp_Pnt myPoint(x,y,z);
2858       BRepClass3d_SolidClassifier scl(shapeToMesh);
2859       scl.Perform(myPoint, 1e-7);
2860       TopAbs_State result = scl.State();
2861       if ( result != TopAbs_IN )
2862         continue;
2863       MESSAGE("Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second);
2864       // X Y Z PHY_SIZE DUMMY_INT
2865       theFile
2866       << x << space
2867       << y << space
2868       << z << space
2869       << vertexIt->second << space
2870       << dummyint << space;
2871       theFile << std::endl;
2872       enfVertexIndexSizeMap[aGhs3dID] = vertexIt->second;
2873       aGhs3dID++;
2874     }
2875   }
2876
2877
2878   std::cout << std::endl;
2879   std::cout << "End writing in 'points' file." << std::endl;
2880
2881   return true;
2882 }
2883
2884 //=======================================================================
2885 //function : readResultFile
2886 //purpose  : readResultFile with geometry
2887 //=======================================================================
2888
2889 static bool readResultFile(const int                       fileOpen,
2890 #ifdef WNT
2891                            const char*                     fileName,
2892 #endif
2893                            GHS3DPlugin_GHS3D*              theAlgo,
2894                            SMESH_MesherHelper&             theHelper,
2895                            TopoDS_Shape                    tabShape[],
2896                            double**                        tabBox,
2897                            const int                       nbShape,
2898                            map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
2899                            std::map <int,int> &            theNodeId2NodeIndexMap,
2900                            bool                            toMeshHoles,
2901                            int                             nbEnforcedVertices,
2902                            int                             nbEnforcedNodes,
2903                            GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges,
2904                            GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles,
2905                            bool                            toMakeGroupsOfDomains)
2906 {
2907   MESSAGE("GHS3DPlugin_GHS3D::readResultFile()");
2908   Kernel_Utils::Localizer loc;
2909   struct stat status;
2910   size_t      length;
2911
2912   std::string tmpStr;
2913
2914   char *ptr, *mapPtr;
2915   char *tetraPtr;
2916   char *shapePtr;
2917
2918   SMESHDS_Mesh* theMeshDS = theHelper.GetMeshDS();
2919
2920   int nbElems, nbNodes, nbInputNodes;
2921   int nbTriangle;
2922   int ID, shapeID, ghs3dShapeID;
2923   int IdShapeRef = 1;
2924   int compoundID =
2925     nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
2926
2927   int *tab, *tabID, *nodeID, *nodeAssigne;
2928   double *coord;
2929   const SMDS_MeshNode **node;
2930
2931   tab    = new int[3];
2932   nodeID = new int[4];
2933   coord  = new double[3];
2934   node   = new const SMDS_MeshNode*[4];
2935
2936   TopoDS_Shape aSolid;
2937   SMDS_MeshNode * aNewNode;
2938   map <int,const SMDS_MeshNode*>::iterator itOnNode;
2939   SMDS_MeshElement* aTet;
2940 #ifdef _DEBUG_
2941   set<int> shapeIDs;
2942 #endif
2943
2944   // Read the file state
2945   fstat(fileOpen, &status);
2946   length   = status.st_size;
2947
2948   // Mapping the result file into memory
2949 #ifdef WNT
2950   HANDLE fd = CreateFile(fileName, GENERIC_READ, FILE_SHARE_READ,
2951                          NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
2952   HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
2953                                         0, (DWORD)length, NULL);
2954   ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
2955 #else
2956   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
2957 #endif
2958   mapPtr = ptr;
2959
2960   ptr      = readMapIntLine(ptr, tab);
2961   tetraPtr = ptr;
2962
2963   nbElems      = tab[0];
2964   nbNodes      = tab[1];
2965   nbInputNodes = tab[2];
2966
2967   nodeAssigne = new int[ nbNodes+1 ];
2968
2969   if (nbShape > 0)
2970     aSolid = tabShape[0];
2971
2972   // Reading the nodeId
2973   for (int i=0; i < 4*nbElems; i++)
2974     strtol(ptr, &ptr, 10);
2975
2976   MESSAGE("nbInputNodes: "<<nbInputNodes);
2977   MESSAGE("nbEnforcedVertices: "<<nbEnforcedVertices);
2978   MESSAGE("nbEnforcedNodes: "<<nbEnforcedNodes);
2979   // Reading the nodeCoor and update the nodeMap
2980   for (int iNode=1; iNode <= nbNodes; iNode++) {
2981     if(theAlgo->computeCanceled())
2982       return false;
2983     for (int iCoor=0; iCoor < 3; iCoor++)
2984       coord[ iCoor ] = strtod(ptr, &ptr);
2985     nodeAssigne[ iNode ] = 1;
2986     if ( iNode > (nbInputNodes-(nbEnforcedVertices+nbEnforcedNodes)) ) {
2987       // Creating SMESH nodes
2988       // - for enforced vertices
2989       // - for vertices of forced edges
2990       // - for ghs3d nodes
2991       nodeAssigne[ iNode ] = 0;
2992       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
2993       theGhs3dIdToNodeMap.insert(theGhs3dIdToNodeMap.end(), make_pair( iNode, aNewNode ));
2994     }
2995   }
2996
2997   // Reading the number of triangles which corresponds to the number of sub-domains
2998   nbTriangle = strtol(ptr, &ptr, 10);
2999
3000   tabID = new int[nbTriangle];
3001   for (int i=0; i < nbTriangle; i++) {
3002     if(theAlgo->computeCanceled())
3003       return false;
3004     tabID[i] = 0;
3005     // find the solid corresponding to GHS3D sub-domain following
3006     // the technique proposed in GHS3D manual in chapter
3007     // "B.4 Subdomain (sub-region) assignment"
3008     int nodeId1 = strtol(ptr, &ptr, 10);
3009     int nodeId2 = strtol(ptr, &ptr, 10);
3010     int nodeId3 = strtol(ptr, &ptr, 10);
3011     if ( nbTriangle > 1 ) {
3012       const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
3013       const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
3014       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
3015       if (!n1 || !n2 || !n3) {
3016         tabID[i] = HOLE_ID;
3017         continue;
3018       }
3019       try {
3020         OCC_CATCH_SIGNALS;
3021 //         tabID[i] = findShapeID( theHelper, n1, n2, n3, toMeshHoles );
3022         tabID[i] = findShapeID( *theHelper.GetMesh(), n1, n2, n3, toMeshHoles );
3023         // -- 0020330: Pb with ghs3d as a submesh
3024         // check that found shape is to be meshed
3025         if ( tabID[i] > 0 ) {
3026           const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] );
3027           bool isToBeMeshed = false;
3028           for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS )
3029             isToBeMeshed = foundShape.IsSame( tabShape[ iS ]);
3030           if ( !isToBeMeshed )
3031             tabID[i] = HOLE_ID;
3032         }
3033         // END -- 0020330: Pb with ghs3d as a submesh
3034 #ifdef _DEBUG_
3035         std::cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << std::endl;
3036 #endif
3037       }
3038       catch ( Standard_Failure & ex)
3039       {
3040 #ifdef _DEBUG_
3041         std::cout << i+1 << " subdomain: Exception caugt: " << ex.GetMessageString() << std::endl;
3042 #endif
3043       }
3044       catch (...) {
3045 #ifdef _DEBUG_
3046         std::cout << i+1 << " subdomain: unknown exception caught " << std::endl;
3047 #endif
3048       }
3049     }
3050   }
3051
3052   shapePtr = ptr;
3053
3054   if ( nbTriangle <= nbShape ) // no holes
3055     toMeshHoles = true; // not avoid creating tetras in holes
3056
3057   // IMP 0022172: [CEA 790] create the groups corresponding to domains
3058   std::vector< std::vector< const SMDS_MeshElement* > > elemsOfDomain( Max( nbTriangle, nbShape ));
3059
3060   // Associating the tetrahedrons to the shapes
3061   shapeID = compoundID;
3062   for (int iElem = 0; iElem < nbElems; iElem++) {
3063     if(theAlgo->computeCanceled())
3064       return false;
3065     for (int iNode = 0; iNode < 4; iNode++) {
3066       ID = strtol(tetraPtr, &tetraPtr, 10);
3067       itOnNode = theGhs3dIdToNodeMap.find(ID);
3068       node[ iNode ] = itOnNode->second;
3069       nodeID[ iNode ] = ID;
3070     }
3071     // We always run GHS3D with "to mesh holes"==TRUE but we must not create
3072     // tetras within holes depending on hypo option,
3073     // so we first check if aTet is inside a hole and then create it
3074     //aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
3075     ghs3dShapeID = 0; // domain ID
3076     if ( nbTriangle > 1 ) {
3077       shapeID = HOLE_ID; // negative shapeID means not to create tetras if !toMeshHoles
3078       ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef;
3079       if ( tabID[ ghs3dShapeID ] == 0 ) {
3080         TopAbs_State state;
3081         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
3082         if ( toMeshHoles || state == TopAbs_IN )
3083           shapeID = theMeshDS->ShapeToIndex( aSolid );
3084         tabID[ ghs3dShapeID ] = shapeID;
3085       }
3086       else
3087         shapeID = tabID[ ghs3dShapeID ];
3088     }
3089     else if ( nbShape > 1 ) {
3090       // Case where nbTriangle == 1 while nbShape == 2 encountered
3091       // with compound of 2 boxes and "To mesh holes"==False,
3092       // so there are no subdomains specified for each tetrahedron.
3093       // Try to guess a solid by a node already bound to shape
3094       shapeID = 0;
3095       for ( int i=0; i<4 && shapeID==0; i++ ) {
3096         if ( nodeAssigne[ nodeID[i] ] == 1 &&
3097              node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
3098              node[i]->getshapeId() > 1 )
3099         {
3100           shapeID = node[i]->getshapeId();
3101         }
3102       }
3103       if ( shapeID==0 ) {
3104         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
3105         shapeID = theMeshDS->ShapeToIndex( aSolid );
3106       }
3107     }
3108     // set new nodes and tetrahedron onto the shape
3109     for ( int i=0; i<4; i++ ) {
3110       if ( nodeAssigne[ nodeID[i] ] == 0 ) {
3111         if ( shapeID != HOLE_ID )
3112           theMeshDS->SetNodeInVolume( node[i], shapeID );
3113         nodeAssigne[ nodeID[i] ] = shapeID;
3114       }
3115     }
3116     if ( toMeshHoles || shapeID != HOLE_ID ) {
3117       aTet = theHelper.AddVolume( node[1], node[0], node[2], node[3],
3118                                   /*id=*/0, /*force3d=*/false);
3119       theMeshDS->SetMeshElementOnShape( aTet, shapeID );
3120       if ( toMakeGroupsOfDomains )
3121       {
3122         if ( int( elemsOfDomain.size() ) < ghs3dShapeID+1 )
3123           elemsOfDomain.resize( ghs3dShapeID+1 );
3124         elemsOfDomain[ ghs3dShapeID ].push_back( aTet );
3125       }
3126     }
3127 #ifdef _DEBUG_
3128     shapeIDs.insert( shapeID );
3129 #endif
3130   }
3131   if ( toMakeGroupsOfDomains )
3132     makeDomainGroups( elemsOfDomain, &theHelper );
3133
3134   // Add enforced elements
3135   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap::const_iterator elemIt;
3136   const SMDS_MeshElement* anElem;
3137   SMDS_ElemIteratorPtr itOnEnfElem;
3138   map<int,int>::const_iterator itOnMap;
3139   shapeID = compoundID;
3140   // Enforced edges
3141   if (theEnforcedEdges.size()) {
3142     (theEnforcedEdges.size() <= 1) ? tmpStr = " enforced edge" : " enforced edges";
3143     std::cout << "Add " << theEnforcedEdges.size() << tmpStr << std::endl;
3144     std::vector< const SMDS_MeshNode* > node( 2 );
3145     // Iterate over the enforced edges
3146     for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
3147       anElem = elemIt->first;
3148       bool addElem = true;
3149       itOnEnfElem = anElem->nodesIterator();
3150       for ( int j = 0; j < 2; ++j ) {
3151         int aNodeID = itOnEnfElem->next()->GetID();
3152         itOnMap = theNodeId2NodeIndexMap.find(aNodeID);
3153         if (itOnMap != theNodeId2NodeIndexMap.end()) {
3154           itOnNode = theGhs3dIdToNodeMap.find((*itOnMap).second);
3155           if (itOnNode != theGhs3dIdToNodeMap.end()) {
3156             node.push_back((*itOnNode).second);
3157 //             shapeID =(*itOnNode).second->getshapeId();
3158           }
3159           else
3160             addElem = false;
3161         }
3162         else
3163           addElem = false;
3164       }
3165       if (addElem) {
3166         aTet = theHelper.AddEdge( node[0], node[1], 0,  false);
3167         theMeshDS->SetMeshElementOnShape( aTet, shapeID );
3168       }
3169     }
3170   }
3171   // Enforced faces
3172   if (theEnforcedTriangles.size()) {
3173     (theEnforcedTriangles.size() <= 1) ? tmpStr = " enforced triangle" : " enforced triangles";
3174     std::cout << "Add " << theEnforcedTriangles.size() << " enforced triangles" << std::endl;
3175     std::vector< const SMDS_MeshNode* > node( 3 );
3176     // Iterate over the enforced triangles
3177     for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
3178       anElem = elemIt->first;
3179       bool addElem = true;
3180       itOnEnfElem = anElem->nodesIterator();
3181       for ( int j = 0; j < 3; ++j ) {
3182         int aNodeID = itOnEnfElem->next()->GetID();
3183         itOnMap = theNodeId2NodeIndexMap.find(aNodeID);
3184         if (itOnMap != theNodeId2NodeIndexMap.end()) {
3185           itOnNode = theGhs3dIdToNodeMap.find((*itOnMap).second);
3186           if (itOnNode != theGhs3dIdToNodeMap.end()) {
3187             node.push_back((*itOnNode).second);
3188 //             shapeID =(*itOnNode).second->getshapeId();
3189           }
3190           else
3191             addElem = false;
3192         }
3193         else
3194           addElem = false;
3195       }
3196       if (addElem) {
3197         aTet = theHelper.AddFace( node[0], node[1], node[2], 0,  false);
3198         theMeshDS->SetMeshElementOnShape( aTet, shapeID );
3199       }
3200     }
3201   }
3202
3203   // Remove nodes of tetras inside holes if !toMeshHoles
3204   if ( !toMeshHoles ) {
3205     itOnNode = theGhs3dIdToNodeMap.find( nbInputNodes );
3206     for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
3207       ID = itOnNode->first;
3208       if ( nodeAssigne[ ID ] == HOLE_ID )
3209         theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
3210     }
3211   }
3212
3213
3214   if ( nbElems ) {
3215     (nbElems <= 1) ? tmpStr = " tetrahedra" : " tetrahedrons";
3216     cout << nbElems << tmpStr << " have been associated to " << nbShape;
3217     (nbShape <= 1) ? tmpStr = " shape" : " shapes";
3218     cout << tmpStr << endl;
3219   }
3220 #ifdef WNT
3221   UnmapViewOfFile(mapPtr);
3222   CloseHandle(hMapObject);
3223   CloseHandle(fd);
3224 #else
3225   munmap(mapPtr, length);
3226 #endif
3227   close(fileOpen);
3228
3229   delete [] tab;
3230   delete [] tabID;
3231   delete [] nodeID;
3232   delete [] coord;
3233   delete [] node;
3234   delete [] nodeAssigne;
3235
3236 #ifdef _DEBUG_
3237   shapeIDs.erase(-1);
3238   if ( shapeIDs.size() != nbShape ) {
3239     (shapeIDs.size() <= 1) ? tmpStr = " solid" : " solids";
3240     std::cout << "Only " << shapeIDs.size() << tmpStr << " of " << nbShape << " found" << std::endl;
3241     for (int i=0; i<nbShape; i++) {
3242       shapeID = theMeshDS->ShapeToIndex( tabShape[i] );
3243       if ( shapeIDs.find( shapeID ) == shapeIDs.end() )
3244         std::cout << "  Solid #" << shapeID << " not found" << std::endl;
3245     }
3246   }
3247 #endif
3248
3249   return true;
3250 }
3251
3252
3253 //=============================================================================
3254 /*!
3255  *Here we are going to use the GHS3D mesher with geometry
3256  */
3257 //=============================================================================
3258
3259 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
3260                                 const TopoDS_Shape& theShape)
3261 {
3262   bool Ok(false);
3263   //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
3264
3265   // we count the number of shapes
3266   // _nbShape = countShape( meshDS, TopAbs_SOLID ); -- 0020330: Pb with ghs3d as a submesh
3267   _nbShape = 0;
3268   TopExp_Explorer expBox ( theShape, TopAbs_SOLID );
3269   for ( ; expBox.More(); expBox.Next() )
3270     _nbShape++;
3271
3272   // create bounding box for every shape inside the compound
3273
3274   int iShape = 0;
3275   TopoDS_Shape* tabShape;
3276   double**      tabBox;
3277   tabShape = new TopoDS_Shape[_nbShape];
3278   tabBox   = new double*[_nbShape];
3279   for (int i=0; i<_nbShape; i++)
3280     tabBox[i] = new double[6];
3281   Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
3282
3283   for (expBox.ReInit(); expBox.More(); expBox.Next()) {
3284     tabShape[iShape] = expBox.Current();
3285     Bnd_Box BoundingBox;
3286     BRepBndLib::Add(expBox.Current(), BoundingBox);
3287     BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
3288     tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
3289     tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
3290     tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
3291     iShape++;
3292   }
3293
3294   // a unique working file name
3295   // to avoid access to the same files by eg different users
3296   _genericName = GHS3DPlugin_Hypothesis::GetFileName(_hyp);
3297   TCollection_AsciiString aGenericName((char*) _genericName.c_str() );
3298
3299   TCollection_AsciiString aResultFileName;
3300   TCollection_AsciiString aLogFileName    = aGenericName + ".log";    // log
3301   // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
3302 //   TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName;
3303 //   TCollection_AsciiString aGenericNameRequired = aGenericName + "_required";
3304 // #ifdef _DEBUG_
3305 //   aGMFFileName    = aGenericName + ".mesh"; // GMF mesh file
3306 //   aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file
3307 //   aRequiredVerticesFileName    = aGenericNameRequired + ".mesh"; // GMF required vertices mesh file
3308 //   aSolFileName    = aGenericName + "_required.sol"; // GMF solution file
3309 // #else
3310 //   aGMFFileName    = aGenericName + ".meshb"; // GMF mesh file
3311 //   aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file
3312 //   aRequiredVerticesFileName    = aGenericNameRequired + ".meshb"; // GMF required vertices mesh file
3313 //   aSolFileName    = aGenericName + "_required.solb"; // GMF solution file
3314 // #endif
3315
3316   TCollection_AsciiString aFacesFileName, aPointsFileName, aBadResFileName, aBbResFileName;
3317
3318   aFacesFileName  = aGenericName + ".faces";  // in faces
3319   aPointsFileName = aGenericName + ".points"; // in points
3320   aResultFileName = aGenericName + ".noboite";// out points and volumes
3321   aBadResFileName = aGenericName + ".boite";  // out bad result
3322   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
3323
3324   // -----------------
3325   // make input files
3326   // -----------------
3327
3328   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
3329   ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
3330
3331   Ok =
3332     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
3333   if (!Ok) {
3334     INFOS( "Can't write into " << aFacesFileName);
3335     return error(SMESH_Comment("Can't write into ") << aFacesFileName);
3336   }
3337
3338   std::map <int,int> aNodeId2NodeIndexMap, aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap;
3339   std::map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
3340   std::map <int, int> nodeID2nodeIndexMap;
3341   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues coordsSizeMap = GHS3DPlugin_Hypothesis::GetEnforcedVerticesCoordsSize(_hyp);
3342   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp);
3343   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp);
3344   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp);
3345 //   TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp);
3346   GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp);
3347
3348   int nbEnforcedVertices = coordsSizeMap.size();
3349   int nbEnforcedNodes = enforcedNodes.size();
3350
3351   std::string tmpStr;
3352   (nbEnforcedNodes <= 1) ? tmpStr = "node" : "nodes";
3353   std::cout << nbEnforcedNodes << " enforced " << tmpStr << " from hypo" << std::endl;
3354   (nbEnforcedVertices <= 1) ? tmpStr = "vertex" : "vertices";
3355   std::cout << nbEnforcedVertices << " enforced " << tmpStr << " from hypo" << std::endl;
3356
3357   SMESH_MesherHelper helper( theMesh );
3358   helper.SetSubShape( theShape );
3359
3360   {
3361     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
3362
3363     // make prisms on quadrangles
3364     if ( theMesh.NbQuadrangles() > 0 )
3365     {
3366       vector<SMESH_ProxyMesh::Ptr> components;
3367       for (expBox.ReInit(); expBox.More(); expBox.Next())
3368       {
3369         if ( _viscousLayersHyp )
3370         {
3371           proxyMesh = _viscousLayersHyp->Compute( theMesh, expBox.Current() );
3372           if ( !proxyMesh )
3373             return false;
3374         }
3375         StdMeshers_QuadToTriaAdaptor* q2t = new StdMeshers_QuadToTriaAdaptor;
3376         q2t->Compute( theMesh, expBox.Current(), proxyMesh.get() );
3377         components.push_back( SMESH_ProxyMesh::Ptr( q2t ));
3378       }
3379       proxyMesh.reset( new SMESH_ProxyMesh( components ));
3380     }
3381     // build viscous layers
3382     else if ( _viscousLayersHyp )
3383     {
3384       proxyMesh = _viscousLayersHyp->Compute( theMesh, theShape );
3385       if ( !proxyMesh )
3386         return false;
3387     }
3388
3389     Ok = (writePoints( aPointsFile, helper,
3390                        aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, aGhs3dIdToNodeMap,
3391                        nodeIDToSizeMap,
3392                        coordsSizeMap, enforcedNodes, enforcedEdges, enforcedTriangles)
3393           &&
3394           writeFaces ( aFacesFile, *proxyMesh, theShape,
3395                        aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap,
3396                        enforcedEdges, enforcedTriangles ));
3397 //     Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(),
3398 //                       helper, *proxyMesh,
3399 //                       aNodeId2NodeIndexMap, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap,
3400 //                       enforcedNodes, enforcedEdges, enforcedTriangles, /*enforcedQuadrangles,*/
3401 //                       coordsSizeMap);
3402   }
3403
3404   // Write aSmdsToGhs3dIdMap to temp file
3405   TCollection_AsciiString aSmdsToGhs3dIdMapFileName;
3406   aSmdsToGhs3dIdMapFileName = aGenericName + ".ids";  // ids relation
3407   ofstream aIdsFile  ( aSmdsToGhs3dIdMapFileName.ToCString()  , ios::out);
3408   Ok = aIdsFile.rdbuf()->is_open();
3409   if (!Ok) {
3410     INFOS( "Can't write into " << aSmdsToGhs3dIdMapFileName);
3411     return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName);
3412   }
3413   INFOS( "Writing ids relation into " << aSmdsToGhs3dIdMapFileName);
3414   aIdsFile << "Smds Ghs3d" << std::endl;
3415   map <int,int>::const_iterator myit;
3416   for (myit=aSmdsToGhs3dIdMap.begin() ; myit != aSmdsToGhs3dIdMap.end() ; ++myit) {
3417     aIdsFile << myit->first << " " << myit->second << std::endl;
3418   }
3419
3420   aIdsFile.close();
3421   aFacesFile.close();
3422   aPointsFile.close();
3423
3424   if ( ! Ok ) {
3425     if ( !_keepFiles ) {
3426 //       removeFile( aGMFFileName );
3427 //       removeFile( aRequiredVerticesFileName );
3428 //       removeFile( aSolFileName );
3429       removeFile( aFacesFileName );
3430       removeFile( aPointsFileName );
3431       removeFile( aSmdsToGhs3dIdMapFileName );
3432     }
3433     return error(COMPERR_BAD_INPUT_MESH);
3434   }
3435   removeFile( aResultFileName ); // needed for boundary recovery module usage
3436
3437   // -----------------
3438   // run ghs3d mesher
3439   // -----------------
3440
3441   TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
3442   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
3443
3444   if ( !_logInStandardOutput )
3445     cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
3446   // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
3447 //   cmd += TCollection_AsciiString(" --in ") + aGenericName;
3448 //   cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired;
3449 //    cmd += TCollection_AsciiString(" --out ") + aResultGMFFileName;
3450 //   cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
3451
3452   std::cout << std::endl;
3453   std::cout << "Ghs3d execution..." << std::endl;
3454   std::cout << cmd << std::endl;
3455
3456   _compute_canceled = false;
3457
3458   system( cmd.ToCString() ); // run
3459
3460   std::cout << std::endl;
3461   std::cout << "End of Ghs3d execution !" << std::endl;
3462
3463   // --------------
3464   // read a result
3465   // --------------
3466
3467   // Mapping the result file
3468
3469   int fileOpen;
3470   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
3471   if ( fileOpen < 0 ) {
3472     std::cout << std::endl;
3473     std::cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << std::endl;
3474     std::cout << "Log: " << aLogFileName << std::endl;
3475     Ok = false;
3476   }
3477   else {
3478     bool toMeshHoles =
3479       _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
3480     const bool toMakeGroupsOfDomains = GHS3DPlugin_Hypothesis::GetToMakeGroupsOfDomains( _hyp );
3481
3482     helper.IsQuadraticSubMesh( theShape );
3483     helper.SetElementsOnShape( false );
3484
3485     Ok = readResultFile( fileOpen,
3486 #ifdef WNT
3487                          aResultFileName.ToCString(),
3488 #endif
3489                          this,
3490                          /*theMesh, */helper, tabShape, tabBox, _nbShape,
3491                          aGhs3dIdToNodeMap, aNodeId2NodeIndexMap,
3492                          toMeshHoles,
3493                          nbEnforcedVertices, nbEnforcedNodes,
3494                          enforcedEdges, enforcedTriangles,
3495                          toMakeGroupsOfDomains );
3496
3497 //       Ok = readGMFFile(
3498 // #ifndef GMF_HAS_SUBDOMAIN_INFO
3499 //                        fileOpen,
3500 // #endif
3501 //                        aGenericName.ToCString(), theMesh,
3502 //                        _nbShape, tabShape, tabBox,
3503 //                        aGhs3dIdToNodeMap, toMeshHoles,
3504 //                        nbEnforcedVertices, nbEnforcedNodes);
3505
3506     removeEmptyGroupsOfDomains( helper.GetMesh(), /*notEmptyAsWell =*/ !toMakeGroupsOfDomains );
3507   }
3508
3509
3510
3511
3512   // ---------------------
3513   // remove working files
3514   // ---------------------
3515
3516   if ( Ok )
3517   {
3518     if ( _removeLogOnSuccess )
3519       removeFile( aLogFileName );
3520
3521     // if ( _hyp && _hyp->GetToMakeGroupsOfDomains() )
3522     //   error( COMPERR_WARNING, "'toMakeGroupsOfDomains' is ignored since the mesh is on shape" );
3523   }
3524   else if ( OSD_File( aLogFileName ).Size() > 0 )
3525   {
3526     // get problem description from the log file
3527     _Ghs2smdsConvertor conv( aGhs3dIdToNodeMap );
3528     storeErrorDescription( aLogFileName, conv );
3529   }
3530   else
3531   {
3532     // the log file is empty
3533     removeFile( aLogFileName );
3534     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
3535     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
3536   }
3537
3538   if ( !_keepFiles ) {
3539     if (! Ok && _compute_canceled)
3540         removeFile( aLogFileName );
3541     removeFile( aFacesFileName );
3542     removeFile( aPointsFileName );
3543     removeFile( aResultFileName );
3544     removeFile( aBadResFileName );
3545     removeFile( aBbResFileName );
3546     removeFile( aSmdsToGhs3dIdMapFileName );
3547   }
3548   std::cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
3549   if ( !Ok )
3550     std::cout << "not ";
3551   std::cout << "treated !" << std::endl;
3552   std::cout << std::endl;
3553
3554   _nbShape = 0;    // re-initializing _nbShape for the next Compute() method call
3555   delete [] tabShape;
3556   delete [] tabBox;
3557
3558   return Ok;
3559 }
3560
3561 //=============================================================================
3562 /*!
3563  *Here we are going to use the GHS3D mesher w/o geometry
3564  */
3565 //=============================================================================
3566 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
3567                                 SMESH_MesherHelper* theHelper)
3568 {
3569   MESSAGE("GHS3DPlugin_GHS3D::Compute()");
3570
3571   //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
3572   TopoDS_Shape theShape = theHelper->GetSubShape();
3573
3574   // a unique working file name
3575   // to avoid access to the same files by eg different users
3576   _genericName = GHS3DPlugin_Hypothesis::GetFileName(_hyp);
3577   TCollection_AsciiString aGenericName((char*) _genericName.c_str() );
3578   TCollection_AsciiString aGenericNameRequired = aGenericName + "_required";
3579
3580   TCollection_AsciiString aLogFileName    = aGenericName + ".log";    // log
3581   TCollection_AsciiString aResultFileName;
3582   bool Ok;
3583
3584   TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName;
3585 //#ifdef _DEBUG_
3586   aGMFFileName              = aGenericName + ".mesh"; // GMF mesh file
3587   aResultFileName           = aGenericName + "Vol.mesh"; // GMF mesh file
3588   aRequiredVerticesFileName = aGenericNameRequired + ".mesh"; // GMF required vertices mesh file
3589   aSolFileName              = aGenericNameRequired + ".sol"; // GMF solution file
3590 //#else
3591 //  aGMFFileName    = aGenericName + ".meshb"; // GMF mesh file
3592 //  aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file
3593 //  aRequiredVerticesFileName    = aGenericNameRequired + ".meshb"; // GMF required vertices mesh file
3594 //  aSolFileName    = aGenericNameRequired + ".solb"; // GMF solution file
3595 //#endif
3596
3597   std::map <int, int> nodeID2nodeIndexMap;
3598   std::map<std::vector<double>, std::string> enfVerticesWithGroup;
3599   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues coordsSizeMap;
3600   TopoDS_Shape GeomShape;
3601 //   TopAbs_ShapeEnum GeomType;
3602   std::vector<double> coords;
3603   gp_Pnt aPnt;
3604   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertex* enfVertex;
3605
3606   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList enfVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
3607   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList::const_iterator enfVerIt = enfVertices.begin();
3608
3609   for ( ; enfVerIt != enfVertices.end() ; ++enfVerIt)
3610   {
3611     enfVertex = (*enfVerIt);
3612 //     if (enfVertex->geomEntry.empty() && enfVertex->coords.size()) {
3613     if (enfVertex->coords.size()) {
3614       coordsSizeMap.insert(make_pair(enfVertex->coords,enfVertex->size));
3615       enfVerticesWithGroup.insert(make_pair(enfVertex->coords,enfVertex->groupName));
3616 //       MESSAGE("enfVerticesWithGroup.insert(make_pair(("<<enfVertex->coords[0]<<","<<enfVertex->coords[1]<<","<<enfVertex->coords[2]<<"),\""<<enfVertex->groupName<<"\"))");
3617     }
3618     else {
3619 //     if (!enfVertex->geomEntry.empty()) {
3620       GeomShape = entryToShape(enfVertex->geomEntry);
3621 //       GeomType = GeomShape.ShapeType();
3622
3623 //       if (!enfVertex->isCompound) {
3624 // //       if (GeomType == TopAbs_VERTEX) {
3625 //         coords.clear();
3626 //         aPnt = BRep_Tool::Pnt(TopoDS::Vertex(GeomShape));
3627 //         coords.push_back(aPnt.X());
3628 //         coords.push_back(aPnt.Y());
3629 //         coords.push_back(aPnt.Z());
3630 //         if (coordsSizeMap.find(coords) == coordsSizeMap.end()) {
3631 //           coordsSizeMap.insert(make_pair(coords,enfVertex->size));
3632 //           enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName));
3633 //         }
3634 //       }
3635 //
3636 //       // Group Management
3637 //       else {
3638 //       if (GeomType == TopAbs_COMPOUND){
3639         for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
3640           coords.clear();
3641           if (it.Value().ShapeType() == TopAbs_VERTEX){
3642             aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value()));
3643             coords.push_back(aPnt.X());
3644             coords.push_back(aPnt.Y());
3645             coords.push_back(aPnt.Z());
3646             if (coordsSizeMap.find(coords) == coordsSizeMap.end()) {
3647               coordsSizeMap.insert(make_pair(coords,enfVertex->size));
3648               enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName));
3649 //               MESSAGE("enfVerticesWithGroup.insert(make_pair(("<<coords[0]<<","<<coords[1]<<","<<coords[2]<<"),\""<<enfVertex->groupName<<"\"))");
3650             }
3651           }
3652         }
3653 //       }
3654     }
3655   }
3656
3657 //   const SMDS_MeshNode* enfNode;
3658   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp);
3659 //   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator enfNodeIt = enforcedNodes.begin();
3660 //   for ( ; enfNodeIt != enforcedNodes.end() ; ++enfNodeIt)
3661 //   {
3662 //     enfNode = enfNodeIt->first;
3663 //     coords.clear();
3664 //     coords.push_back(enfNode->X());
3665 //     coords.push_back(enfNode->Y());
3666 //     coords.push_back(enfNode->Z());
3667 //     if (enfVerticesWithGro
3668 //       enfVerticesWithGroup.insert(make_pair(coords,enfNodeIt->second));
3669 //   }
3670
3671
3672   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp);
3673   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp);
3674 //   TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp);
3675   GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp);
3676
3677   std::string tmpStr;
3678
3679   int nbEnforcedVertices = coordsSizeMap.size();
3680   int nbEnforcedNodes = enforcedNodes.size();
3681   (nbEnforcedNodes <= 1) ? tmpStr = "node" : tmpStr = "nodes";
3682   std::cout << nbEnforcedNodes << " enforced " << tmpStr << " from hypo" << std::endl;
3683   (nbEnforcedVertices <= 1) ? tmpStr = "vertex" : tmpStr = "vertices";
3684   std::cout << nbEnforcedVertices << " enforced " << tmpStr << " from hypo" << std::endl;
3685
3686   std::vector <const SMDS_MeshNode*> aNodeByGhs3dId, anEnforcedNodeByGhs3dId;
3687   std::map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
3688   std::vector<std::string> aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId;
3689   {
3690     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
3691     if ( theMesh.NbQuadrangles() > 0 )
3692     {
3693       StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor;
3694       aQuad2Trias->Compute( theMesh );
3695       proxyMesh.reset( aQuad2Trias );
3696     }
3697
3698     Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(),
3699                       *proxyMesh, &theMesh,
3700                       aNodeByGhs3dId, aNodeToGhs3dIdMap,
3701                       aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId,
3702                       enforcedNodes, enforcedEdges, enforcedTriangles,
3703                       enfVerticesWithGroup, coordsSizeMap);
3704   }
3705
3706   // -----------------
3707   // run ghs3d mesher
3708   // -----------------
3709
3710   TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str());
3711
3712   cmd += TCollection_AsciiString(" --in ") + aGMFFileName;
3713   if ( nbEnforcedVertices + nbEnforcedNodes)
3714     cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired;
3715   cmd += TCollection_AsciiString(" --out ") + aResultFileName;
3716   if ( !_logInStandardOutput )
3717     cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
3718
3719   std::cout << std::endl;
3720   std::cout << "Ghs3d execution..." << std::endl;
3721   std::cout << cmd << std::endl;
3722
3723   _compute_canceled = false;
3724
3725   system( cmd.ToCString() ); // run
3726
3727   std::cout << std::endl;
3728   std::cout << "End of Ghs3d execution !" << std::endl;
3729
3730   // --------------
3731   // read a result
3732   // --------------
3733   GHS3DPlugin_Hypothesis::TSetStrings groupsToRemove = GHS3DPlugin_Hypothesis::GetGroupsToRemove(_hyp);
3734   const bool toMakeGroupsOfDomains = GHS3DPlugin_Hypothesis::GetToMakeGroupsOfDomains( _hyp );
3735
3736   Ok = readGMFFile(aResultFileName.ToCString(),
3737                    this,
3738                    theHelper, theShape, aNodeByGhs3dId, aNodeToGhs3dIdMap,
3739                    aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId,
3740                    groupsToRemove, toMakeGroupsOfDomains);
3741
3742   updateMeshGroups(theHelper->GetMesh(), groupsToRemove);
3743   removeEmptyGroupsOfDomains( theHelper->GetMesh(), /*notEmptyAsWell =*/ !toMakeGroupsOfDomains );
3744
3745   if ( Ok ) {
3746     GHS3DPlugin_Hypothesis* that = (GHS3DPlugin_Hypothesis*)this->_hyp;
3747     if (that)
3748       that->ClearGroupsToRemove();
3749   }
3750   // ---------------------
3751   // remove working files
3752   // ---------------------
3753
3754   if ( Ok )
3755   {
3756     if ( _removeLogOnSuccess )
3757       removeFile( aLogFileName );
3758
3759     //if ( !toMakeGroupsOfDomains && _hyp && _hyp->GetToMakeGroupsOfDomains() )
3760     //error( COMPERR_WARNING, "'toMakeGroupsOfDomains' is ignored since 'toMeshHoles' is OFF." );
3761   }
3762   else if ( OSD_File( aLogFileName ).Size() > 0 )
3763   {
3764     // get problem description from the log file
3765     _Ghs2smdsConvertor conv( aNodeByGhs3dId );
3766     storeErrorDescription( aLogFileName, conv );
3767   }
3768   else {
3769     // the log file is empty
3770     removeFile( aLogFileName );
3771     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
3772     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
3773   }
3774
3775   if ( !_keepFiles )
3776   {
3777     if (! Ok && _compute_canceled)
3778         removeFile( aLogFileName );
3779     removeFile( aGMFFileName );
3780     removeFile( aResultFileName );
3781     removeFile( aRequiredVerticesFileName );
3782     removeFile( aSolFileName );
3783   }
3784   return Ok;
3785 }
3786
3787 void GHS3DPlugin_GHS3D::CancelCompute()
3788 {
3789   _compute_canceled = true;
3790 #ifdef WNT
3791 #else
3792   std::string cmd = "ps xo pid,args | grep " + _genericName;
3793   //cmd += " | grep -e \"^ *[0-9]\\+ \\+" + GHS3DPlugin_Hypothesis::GetExeName() + "\"";
3794   cmd += " | awk '{print $1}' | xargs kill -9 > /dev/null 2>&1";
3795   system( cmd.c_str() );
3796 #endif
3797 }
3798
3799 //================================================================================
3800 /*!
3801  * \brief Provide human readable text by error code reported by ghs3d
3802  */
3803 //================================================================================
3804
3805 static const char* translateError(const int errNum)
3806 {
3807   switch ( errNum ) {
3808   case 0:
3809     return "The surface mesh includes a face of type other than edge, "
3810       "triangle or quadrilateral. This face type is not supported.";
3811   case 1:
3812     return "Not enough memory for the face table.";
3813   case 2:
3814     return "Not enough memory.";
3815   case 3:
3816     return "Not enough memory.";
3817   case 4:
3818     return "Face is ignored.";
3819   case 5:
3820     return "End of file. Some data are missing in the file.";
3821   case 6:
3822     return "Read error on the file. There are wrong data in the file.";
3823   case 7:
3824     return "the metric file is inadequate (dimension other than 3).";
3825   case 8:
3826     return "the metric file is inadequate (values not per vertices).";
3827   case 9:
3828     return "the metric file contains more than one field.";
3829   case 10:
3830     return "the number of values in the \".bb\" (metric file) is incompatible with the expected"
3831       "value of number of mesh vertices in the \".noboite\" file.";
3832   case 12:
3833     return "Too many sub-domains.";
3834   case 13:
3835     return "the number of vertices is negative or null.";
3836   case 14:
3837     return "the number of faces is negative or null.";
3838   case 15:
3839     return "A face has a null vertex.";
3840   case 22:
3841     return "incompatible data.";
3842   case 131:
3843     return "the number of vertices is negative or null.";
3844   case 132:
3845     return "the number of vertices is negative or null (in the \".mesh\" file).";
3846   case 133:
3847     return "the number of faces is negative or null.";
3848   case 1000:
3849     return "A face appears more than once in the input surface mesh.";
3850   case 1001:
3851     return "An edge appears more than once in the input surface mesh.";
3852   case 1002:
3853     return "A face has a vertex negative or null.";
3854   case 1003:
3855     return "NOT ENOUGH MEMORY.";
3856   case 2000:
3857     return "Not enough available memory.";
3858   case 2002:
3859     return "Some initial points cannot be inserted. The surface mesh is probably very bad "
3860       "in terms of quality or the input list of points is wrong.";
3861   case 2003:
3862     return "Some vertices are too close to one another or coincident.";
3863   case 2004:
3864     return "Some vertices are too close to one another or coincident.";
3865   case 2012:
3866     return "A vertex cannot be inserted.";
3867   case 2014:
3868     return "There are at least two points considered as coincident.";
3869   case 2103:
3870     return "Some vertices are too close to one another or coincident.";
3871   case 3000:
3872     return "The surface mesh regeneration step has failed.";
3873   case 3009:
3874     return "Constrained edge cannot be enforced.";
3875   case 3019:
3876     return "Constrained face cannot be enforced.";
3877   case 3029:
3878     return "Missing faces.";
3879   case 3100:
3880     return "No guess to start the definition of the connected component(s).";
3881   case 3101:
3882     return "The surface mesh includes at least one hole. The domain is not well defined.";
3883   case 3102:
3884     return "Impossible to define a component.";
3885   case 3103:
3886     return "The surface edge intersects another surface edge.";
3887   case 3104:
3888     return "The surface edge intersects the surface face.";
3889   case 3105:
3890     return "One boundary point lies within a surface face.";
3891   case 3106:
3892     return "One surface edge intersects a surface face.";
3893   case 3107:
3894     return "One boundary point lies within a surface edge.";
3895   case 3108:
3896     return "Insufficient memory ressources detected due to a bad quality surface mesh leading "
3897       "to too many swaps.";
3898   case 3109:
3899     return "Edge is unique (i.e., bounds a hole in the surface).";
3900   case 3122:
3901     return "Presumably, the surface mesh is not compatible with the domain being processed.";
3902   case 3123:
3903     return "Too many components, too many sub-domain.";
3904   case 3209:
3905     return "The surface mesh includes at least one hole. "
3906       "Therefore there is no domain properly defined.";
3907   case 3300:
3908     return "Statistics.";
3909   case 3400:
3910     return "Statistics.";
3911   case 3500:
3912     return "Warning, it is dramatically tedious to enforce the boundary items.";
3913   case 4000:
3914     return "Not enough memory at this time, nevertheless, the program continues. "
3915       "The expected mesh will be correct but not really as large as required.";
3916   case 4002:
3917     return "see above error code, resulting quality may be poor.";
3918   case 4003:
3919     return "Not enough memory at this time, nevertheless, the program continues (warning).";
3920   case 8000:
3921     return "Unknown face type.";
3922   case 8005:
3923   case 8006:
3924     return "End of file. Some data are missing in the file.";
3925   case 9000:
3926     return "A too small volume element is detected.";
3927   case 9001:
3928     return "There exists at least a null or negative volume element.";
3929   case 9002:
3930     return "There exist null or negative volume elements.";
3931   case 9003:
3932     return "A too small volume element is detected. A face is considered being degenerated.";
3933   case 9100:
3934     return "Some element is suspected to be very bad shaped or wrong.";
3935   case 9102:
3936     return "A too bad quality face is detected. This face is considered degenerated.";
3937   case 9112:
3938     return "A too bad quality face is detected. This face is degenerated.";
3939   case 9122:
3940     return "Presumably, the surface mesh is not compatible with the domain being processed.";
3941   case 9999:
3942     return "Abnormal error occured, contact hotline.";
3943   case 23600:
3944     return "Not enough memory for the face table.";
3945   case 23601:
3946     return "The algorithm cannot run further. "
3947       "The surface mesh is probably very bad in terms of quality.";
3948   case 23602:
3949     return "Bad vertex number.";
3950   case 1001200:
3951     return "Cannot close mesh file NomFil.";
3952   case 1002010:
3953     return "There are wrong data.";
3954   case 1002120:
3955     return "The number of faces is negative or null.";
3956   case 1002170:
3957     return "The number of vertices is negative or null in the '.sol' file.";
3958   case 1002190:
3959     return "The number of tetrahedra is negative or null.";
3960   case 1002210:
3961     return "The number of vertices is negative or null.";
3962   case 1002211:
3963     return "A face has a vertex negative or null.";
3964   case 1002270:
3965     return "The field is not a size in file NomFil.";
3966   case 1002280:
3967     return "A count is wrong in the enclosing box in the .boite.mesh input "
3968       "file (option '--read_boite').";
3969   case 1002290:
3970     return "A tetrahedron has a vertex with a negative number.";
3971   case 1002300:
3972     return "the 'MeshVersionFormatted' is not 1 or 2 in the '.mesh' file or the '.sol'.";
3973  case 1002370:
3974    return "The number of values in the '.sol' (metric file) is incompatible with "
3975      "the expected value of number of mesh vertices in the '.mesh' file.";
3976   case 1003000:
3977     return "Not enough memory.";
3978   case 1003020:
3979     return "Not enough memory for the face table.";
3980   case 1003050:
3981     return "Insufficient memory ressources detected due to a bad quality "
3982       "surface mesh leading to too many swaps.";
3983   case 1005010:
3984     return "The surface coordinates of a vertex are differing from the "
3985       "volume coordinates, probably due to a precision problem.";
3986   case 1005050:
3987     return "Invalid dimension. Dimension 3 expected.";
3988   case 1005100:
3989     return "A point has a tag 0. This point is probably outside the domain which has been meshed.";
3990   case 1005103:
3991     return "The vertices of an element are too close to one another or coincident.";
3992   case 1005104:
3993     return "There are at least two points whose distance is very small, and considered as coincident.";
3994   case 1005105:
3995     return "Two vertices are too close to one another or coincident.";
3996   case 1005106:
3997     return "A vertex cannot be inserted.";
3998   case 1005107:
3999     return "Two vertices are too close to one another or coincident. Note : When "
4000       "this error occurs during the overconstrained processing phase, this is only "
4001       "a warning which means that it is difficult to break some overconstrained facets.";
4002   case 1005110:
4003     return "Two surface edges are intersecting.";
4004   case 1005120:
4005     return "A surface edge intersects a surface face.";
4006   case 1005150:
4007     return "A boundary point lies within a surface face.";
4008   case 1005160:
4009     return "A boundary point lies within a surface edge.";
4010   case 1005200:
4011     return "A surface mesh appears more than once in the input surface mesh.";
4012   case 1005210:
4013     return "An edge appears more than once in the input surface mesh.";
4014   case 1005225:
4015     return "Surface with unvalid triangles.";
4016   case 1005270:
4017     return "The metric in the '.sol' file contains more than one field.";
4018   case 1005300:
4019     return "The surface mesh includes at least one hole. The domain is not well defined.";
4020   case 1005301:
4021     return "Presumably, the surface mesh is not compatible with the domain being processed (warning).";
4022   case 1005302:
4023     return "Probable faces overlapping somewher.";
4024   case 1005320:
4025     return "The quadratic version does not work with prescribed free edges.";
4026   case 1005321:
4027     return "The quadratic version does not work with a volume mesh.";
4028   case 1005370:
4029     return "The metric in the '.sol' file is inadequate (values not per vertices).";
4030   case 1005371:
4031     return "The number of vertices in the '.sol' is different from the one in the "
4032       "'.mesh' file for the required vertices (option '--required_vertices').";
4033   case 1005372:
4034     return "More than one type in file NomFil. The type must be equal to 1 in the '.sol'"
4035       "for the required vertices (option '--required_vertices').";
4036   case 1005515:
4037     return "Bad vertex number.";
4038   case 1005560:
4039     return "No guess to start the definition of the connected component(s).";
4040   case 1005602:
4041     return "Some initial points cannot be inserted.";
4042   case 1005620:
4043     return "A too bad quality face is detected. This face is considered degenerated.";
4044   case 1005621:
4045     return "A too bad quality face is detected. This face is degenerated.";
4046   case 1005622:
4047     return "The algorithm cannot run further.";
4048   case 1005690:
4049     return "A too small volume element is detected.";
4050   case 1005691:
4051     return "A tetrahedra is suspected to be very bad shaped or wrong.";
4052   case 1005692:
4053     return "There is at least a null or negative volume element. The resulting mesh"
4054       "may be inappropriate.";
4055   case 1005693:
4056     return "There are some null or negative volume element. The resulting mesh may"
4057       "be inappropriate.";
4058   case 1005820:
4059     return "An edge is unique (i.e., bounds a hole in the surface).";
4060   case 1007000:
4061     return "Abnormal or internal error.";
4062   case 1007010:
4063     return "Too many components with respect to too many sub-domain.";
4064   case 1007400:
4065     return "An internal error has been encountered or a signal has been received. "
4066       "Current mesh will not be saved.";
4067   case 1008491:
4068     return "Impossible to define a component.";
4069   case 1008410:
4070     return "There are some overconstrained edges.";
4071   case 1008420:
4072     return "There are some overconstrained facets.";
4073   case 1008422:
4074     return "Give the number of missing faces (information given when regeneration phase failed).";
4075   case 1008423:
4076     return "A constrained face cannot be enforced (information given when regeneration phase failed).";
4077   case 1008441:
4078     return "A constrained edge cannot be enforced.";
4079   case 1008460:
4080     return "It is dramatically tedious to enforce the boundary items.";
4081   case 1008480:
4082     return "The surface mesh regeneration step has failed. A .boite.mesh and .boite.map files are created.";
4083   case 1008490:
4084     return "Invalid resulting mesh.";
4085   case 1008495:
4086     return "P2 correction not successful.";
4087   case 1009000:
4088     return "Program has received an interruption or a termination signal sent by the "
4089       "user or the system administrator. Current mesh will not be saved.";
4090   }
4091   return "";
4092 }
4093
4094 //================================================================================
4095 /*!
4096  * \brief Retrieve from a string given number of integers
4097  */
4098 //================================================================================
4099
4100 static char* getIds( char* ptr, int nbIds, vector<int>& ids )
4101 {
4102   ids.clear();
4103   ids.reserve( nbIds );
4104   while ( nbIds )
4105   {
4106     while ( !isdigit( *ptr )) ++ptr;
4107     if ( ptr[-1] == '-' ) --ptr;
4108     ids.push_back( strtol( ptr, &ptr, 10 ));
4109     --nbIds;
4110   }
4111   return ptr;
4112 }
4113
4114 //================================================================================
4115 /*!
4116  * \brief Retrieve problem description form a log file
4117  *  \retval bool - always false
4118  */
4119 //================================================================================
4120
4121 bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& logFile,
4122                                               const _Ghs2smdsConvertor &     toSmdsConvertor )
4123 {
4124   if(_compute_canceled)
4125     return error(SMESH_Comment("interruption initiated by user"));
4126   // open file
4127 #ifdef WNT
4128   int file = ::_open (logFile.ToCString(), _O_RDONLY|_O_BINARY);
4129 #else
4130   int file = ::open (logFile.ToCString(), O_RDONLY);
4131 #endif
4132   if ( file < 0 )
4133     return error( SMESH_Comment("See ") << logFile << " for problem description");
4134
4135   // get file size
4136   off_t length = lseek( file, 0, SEEK_END);
4137   lseek( file, 0, SEEK_SET);
4138
4139   // read file
4140   vector< char > buf( length );
4141   int nBytesRead = ::read (file, & buf[0], length);
4142   ::close (file);
4143   char* ptr = & buf[0];
4144   char* bufEnd = ptr + nBytesRead;
4145
4146   SMESH_Comment errDescription;
4147
4148   enum { NODE = 1, EDGE, TRIA, VOL, SKIP_ID = 1 };
4149
4150   // look for MeshGems version
4151   // Since "MG-TETRA -- MeshGems 1.1-3 (January, 2013)" error codes change.
4152   // To discriminate old codes from new ones we add 1000000 to the new codes.
4153   // This way value of the new codes is same as absolute value of codes printed
4154   // in the log after "MGMESSAGE" string.
4155   int versionAddition = 0;
4156   {
4157     char* verPtr = ptr;
4158     while ( ++verPtr < bufEnd )
4159     {
4160       if ( strncmp( verPtr, "MG-TETRA -- MeshGems ", 21 ) != 0 )
4161         continue;
4162       if ( strcmp( verPtr, "MG-TETRA -- MeshGems 1.1-3 " ) >= 0 )
4163         versionAddition = 1000000;
4164       ptr = verPtr;
4165       break;
4166     }
4167   }
4168
4169   // look for errors "ERR #"
4170
4171   set<string> foundErrorStr; // to avoid reporting same error several times
4172   set<int>    elemErrorNums; // not to report different types of errors with bad elements
4173   while ( ++ptr < bufEnd )
4174   {
4175     if ( strncmp( ptr, "ERR ", 4 ) != 0 )
4176       continue;
4177
4178     list<const SMDS_MeshElement*> badElems;
4179     vector<int> nodeIds;
4180
4181     ptr += 4;
4182     char* errBeg = ptr;
4183     int   errNum = strtol(ptr, &ptr, 10) + versionAddition;
4184     // we treat errors enumerated in [SALOME platform 0019316] issue
4185     // and all errors from a new (Release 1.1) MeshGems User Manual
4186     switch ( errNum ) {
4187     case 0015: // The face number (numfac) with vertices (f 1, f 2, f 3) has a null vertex.
4188     case 1005620 : // a too bad quality face is detected. This face is considered degenerated.
4189       ptr = getIds(ptr, SKIP_ID, nodeIds);
4190       ptr = getIds(ptr, TRIA, nodeIds);
4191       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4192       break;
4193     case 1005621 : // a too bad quality face is detected. This face is degenerated.
4194       // hence the is degenerated it is invisible, add its edges in addition
4195       ptr = getIds(ptr, SKIP_ID, nodeIds);
4196       ptr = getIds(ptr, TRIA, nodeIds);
4197       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4198       {
4199         vector<int> edgeNodes( nodeIds.begin(), --nodeIds.end() ); // 01
4200         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
4201         edgeNodes[1] = nodeIds[2]; // 02
4202         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
4203         edgeNodes[0] = nodeIds[1]; // 12
4204       }
4205       break;
4206     case 1000: // Face (f 1, f 2, f 3) appears more than once in the input surface mesh.
4207       // ERR  1000 :  1 3 2
4208     case 1002: // Face (f 1, f 2, f 3) has a vertex negative or null
4209     case 3019: // Constrained face (f 1, f 2, f 3) cannot be enforced
4210     case 1002211: // a face has a vertex negative or null.
4211     case 1005200 : // a surface mesh appears more than once in the input surface mesh.
4212     case 1008423 : // a constrained face cannot be enforced (regeneration phase failed).
4213       ptr = getIds(ptr, TRIA, nodeIds);
4214       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4215       break;
4216     case 1001: // Edge (e1, e2) appears more than once in the input surface mesh
4217     case 3009: // Constrained edge (e1, e2) cannot be enforced (warning).
4218       // ERR  3109 :  EDGE  5 6 UNIQUE
4219     case 3109: // Edge (e1, e2) is unique (i.e., bounds a hole in the surface)
4220     case 1005210 : // an edge appears more than once in the input surface mesh.
4221     case 1005820 : // an edge is unique (i.e., bounds a hole in the surface).
4222     case 1008441 : // a constrained edge cannot be enforced.
4223       ptr = getIds(ptr, EDGE, nodeIds);
4224       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4225       break;
4226     case 2004: // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
4227     case 2014: // at least two points whose distance is dist, i.e., considered as coincident
4228     case 2103: // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
4229       // ERR  2103 :  16 WITH  3
4230     case 1005105 : // two vertices are too close to one another or coincident.
4231     case 1005107: // Two vertices are too close to one another or coincident.
4232       ptr = getIds(ptr, NODE, nodeIds);
4233       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4234       ptr = getIds(ptr, NODE, nodeIds);
4235       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4236       break;
4237     case 2012: // Vertex v1 cannot be inserted (warning).
4238     case 1005106 : // a vertex cannot be inserted.
4239       ptr = getIds(ptr, NODE, nodeIds);
4240       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4241       break;
4242     case 3103: // The surface edge (e1, e2) intersects another surface edge (e3, e4)
4243     case 1005110 : // two surface edges are intersecting.
4244       // ERR  3103 :  1 2 WITH  7 3
4245       ptr = getIds(ptr, EDGE, nodeIds);
4246       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4247       ptr = getIds(ptr, EDGE, nodeIds);
4248       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4249       break;
4250     case 3104: // The surface edge (e1, e2) intersects the surface face (f 1, f 2, f 3)
4251       // ERR  3104 :  9 10 WITH  1 2 3
4252     case 3106: // One surface edge (say e1, e2) intersects a surface face (f 1, f 2, f 3)
4253     case 1005120 : // a surface edge intersects a surface face.
4254       ptr = getIds(ptr, EDGE, nodeIds);
4255       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4256       ptr = getIds(ptr, TRIA, nodeIds);
4257       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4258       break;
4259     case 3105: // One boundary point (say p1) lies within a surface face (f 1, f 2, f 3)
4260       // ERR  3105 :  8 IN  2 3 5
4261     case 1005150 : // a boundary point lies within a surface face.
4262       ptr = getIds(ptr, NODE, nodeIds);
4263       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4264       ptr = getIds(ptr, TRIA, nodeIds);
4265       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4266       break;
4267     case 3107: // One boundary point (say p1) lies within a surface edge (e1, e2) (stop).
4268       // ERR  3107 :  2 IN  4 1
4269     case 1005160 : // a boundary point lies within a surface edge.
4270       ptr = getIds(ptr, NODE, nodeIds);
4271       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4272       ptr = getIds(ptr, EDGE, nodeIds);
4273       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4274       break;
4275     case 9000: // ERR  9000
4276       //  ELEMENT  261 WITH VERTICES :  7 396 -8 242
4277       //  VOLUME   : -1.11325045E+11 W.R.T. EPSILON   0.
4278       // A too small volume element is detected. Are reported the index of the element,
4279       // its four vertex indices, its volume and the tolerance threshold value
4280       ptr = getIds(ptr, SKIP_ID, nodeIds);
4281       ptr = getIds(ptr, VOL, nodeIds);
4282       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4283       // even if all nodes found, volume it most probably invisible,
4284       // add its faces to demonstrate it anyhow
4285       {
4286         vector<int> faceNodes( nodeIds.begin(), --nodeIds.end() ); // 012
4287         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
4288         faceNodes[2] = nodeIds[3]; // 013
4289         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
4290         faceNodes[1] = nodeIds[2]; // 023
4291         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
4292         faceNodes[0] = nodeIds[1]; // 123
4293         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
4294       }
4295       break;
4296     case 9001: // ERR  9001
4297       //  %% NUMBER OF NEGATIVE VOLUME TETS  :  1
4298       //  %% THE LARGEST NEGATIVE TET        :   1.75376581E+11
4299       //  %%  NUMBER OF NULL VOLUME TETS     :  0
4300       // There exists at least a null or negative volume element
4301       break;
4302     case 9002:
4303       // There exist n null or negative volume elements
4304       break;
4305     case 9003:
4306       // A too small volume element is detected
4307       break;
4308     case 9102:
4309       // A too bad quality face is detected. This face is considered degenerated,
4310       // its index, its three vertex indices together with its quality value are reported
4311       break; // same as next
4312     case 9112: // ERR  9112
4313       //  FACE   2 WITH VERTICES :  4 2 5
4314       //  SMALL INRADIUS :   0.
4315       // A too bad quality face is detected. This face is degenerated,
4316       // its index, its three vertex indices together with its inradius are reported
4317       ptr = getIds(ptr, SKIP_ID, nodeIds);
4318       ptr = getIds(ptr, TRIA, nodeIds);
4319       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4320       // add triangle edges as it most probably has zero area and hence invisible
4321       {
4322         vector<int> edgeNodes(2);
4323         edgeNodes[0] = nodeIds[0]; edgeNodes[1] = nodeIds[1]; // 0-1
4324         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
4325         edgeNodes[1] = nodeIds[2]; // 0-2
4326         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
4327         edgeNodes[0] = nodeIds[1]; // 1-2
4328         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
4329       }
4330       break;
4331     case 1005103 : // the vertices of an element are too close to one another or coincident.
4332       ptr = getIds(ptr, TRIA, nodeIds);
4333       if ( nodeIds.back() == 0 ) // index of the third vertex of the element (0 for an edge)
4334         nodeIds.resize( EDGE );
4335       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
4336       break;
4337     }
4338
4339     bool isNewError = foundErrorStr.insert( string( errBeg, ptr )).second;
4340     if ( !isNewError )
4341       continue; // not to report same error several times
4342
4343 //     const SMDS_MeshElement* nullElem = 0;
4344 //     bool allElemsOk = ( find( badElems.begin(), badElems.end(), nullElem) == badElems.end());
4345
4346 //     if ( allElemsOk && !badElems.empty() && !elemErrorNums.empty() ) {
4347 //       bool oneMoreErrorType = elemErrorNums.insert( errNum ).second;
4348 //       if ( oneMoreErrorType )
4349 //         continue; // not to report different types of errors with bad elements
4350 //     }
4351
4352     // store bad elements
4353     //if ( allElemsOk ) {
4354       list<const SMDS_MeshElement*>::iterator elem = badElems.begin();
4355       for ( ; elem != badElems.end(); ++elem )
4356         addBadInputElement( *elem );
4357       //}
4358
4359     // make error text
4360     string text = translateError( errNum );
4361     if ( errDescription.find( text ) == text.npos ) {
4362       if ( !errDescription.empty() )
4363         errDescription << "\n";
4364       errDescription << text;
4365     }
4366
4367   } // end while
4368
4369   if ( errDescription.empty() ) { // no errors found
4370     char msgLic1[] = "connection to server failed";
4371     char msgLic2[] = " Dlim ";
4372     if ( search( &buf[0], bufEnd, msgLic1, msgLic1 + strlen(msgLic1)) != bufEnd ||
4373          search( &buf[0], bufEnd, msgLic2, msgLic2 + strlen(msgLic2)) != bufEnd )
4374       errDescription << "Licence problems.";
4375     else
4376     {
4377       char msg2[] = "SEGMENTATION FAULT";
4378       if ( search( &buf[0], bufEnd, msg2, msg2 + strlen(msg2)) != bufEnd )
4379         errDescription << "ghs3d: SEGMENTATION FAULT. ";
4380     }
4381   }
4382
4383   if ( errDescription.empty() )
4384     errDescription << "See " << logFile << " for problem description";
4385   else
4386     errDescription << "\nSee " << logFile << " for more information";
4387
4388   return error( errDescription );
4389 }
4390
4391 //================================================================================
4392 /*!
4393  * \brief Creates _Ghs2smdsConvertor
4394  */
4395 //================================================================================
4396
4397 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const map <int,const SMDS_MeshNode*> & ghs2NodeMap)
4398   :_ghs2NodeMap( & ghs2NodeMap ), _nodeByGhsId( 0 )
4399 {
4400 }
4401
4402 //================================================================================
4403 /*!
4404  * \brief Creates _Ghs2smdsConvertor
4405  */
4406 //================================================================================
4407
4408 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const vector <const SMDS_MeshNode*> &  nodeByGhsId)
4409   : _ghs2NodeMap( 0 ), _nodeByGhsId( &nodeByGhsId )
4410 {
4411 }
4412
4413 //================================================================================
4414 /*!
4415  * \brief Return SMDS element by ids of GHS3D nodes
4416  */
4417 //================================================================================
4418
4419 const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector<int>& ghsNodes) const
4420 {
4421   size_t nbNodes = ghsNodes.size();
4422   vector<const SMDS_MeshNode*> nodes( nbNodes, 0 );
4423   for ( size_t i = 0; i < nbNodes; ++i ) {
4424     int ghsNode = ghsNodes[ i ];
4425     if ( _ghs2NodeMap ) {
4426       map <int,const SMDS_MeshNode*>::const_iterator in = _ghs2NodeMap->find( ghsNode);
4427       if ( in == _ghs2NodeMap->end() )
4428         return 0;
4429       nodes[ i ] = in->second;
4430     }
4431     else {
4432       if ( ghsNode < 1 || ghsNode > _nodeByGhsId->size() )
4433         return 0;
4434       nodes[ i ] = (*_nodeByGhsId)[ ghsNode-1 ];
4435     }
4436   }
4437   if ( nbNodes == 1 )
4438     return nodes[0];
4439
4440   if ( nbNodes == 2 ) {
4441     const SMDS_MeshElement* edge= SMDS_Mesh::FindEdge( nodes[0], nodes[1] );
4442     if ( !edge )
4443       edge = new SMDS_LinearEdge( nodes[0], nodes[1] );
4444     return edge;
4445   }
4446   if ( nbNodes == 3 ) {
4447     const SMDS_MeshElement* face = SMDS_Mesh::FindFace( nodes );
4448     if ( !face )
4449       face = new SMDS_FaceOfNodes( nodes[0], nodes[1], nodes[2] );
4450     return face;
4451   }
4452   if ( nbNodes == 4 )
4453     return new SMDS_VolumeOfNodes( nodes[0], nodes[1], nodes[2], nodes[3] );
4454
4455   return 0;
4456 }
4457
4458
4459 //=============================================================================
4460 /*!
4461  *
4462  */
4463 //=============================================================================
4464 bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh,
4465                                  const TopoDS_Shape& aShape,
4466                                  MapShapeNbElems& aResMap)
4467 {
4468   int nbtri = 0, nbqua = 0;
4469   double fullArea = 0.0;
4470   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
4471     TopoDS_Face F = TopoDS::Face( exp.Current() );
4472     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
4473     MapShapeNbElemsItr anIt = aResMap.find(sm);
4474     if( anIt==aResMap.end() ) {
4475       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
4476       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
4477                                             "Submesh can not be evaluated",this));
4478       return false;
4479     }
4480     std::vector<int> aVec = (*anIt).second;
4481     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
4482     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
4483     GProp_GProps G;
4484     BRepGProp::SurfaceProperties(F,G);
4485     double anArea = G.Mass();
4486     fullArea += anArea;
4487   }
4488
4489   // collect info from edges
4490   int nb0d_e = 0, nb1d_e = 0;
4491   bool IsQuadratic = false;
4492   bool IsFirst = true;
4493   TopTools_MapOfShape tmpMap;
4494   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
4495     TopoDS_Edge E = TopoDS::Edge(exp.Current());
4496     if( tmpMap.Contains(E) )
4497       continue;
4498     tmpMap.Add(E);
4499     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
4500     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
4501     std::vector<int> aVec = (*anIt).second;
4502     nb0d_e += aVec[SMDSEntity_Node];
4503     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
4504     if(IsFirst) {
4505       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
4506       IsFirst = false;
4507     }
4508   }
4509   tmpMap.Clear();
4510
4511   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
4512
4513   GProp_GProps G;
4514   BRepGProp::VolumeProperties(aShape,G);
4515   double aVolume = G.Mass();
4516   double tetrVol = 0.1179*ELen*ELen*ELen;
4517   double CoeffQuality = 0.9;
4518   int nbVols = int(aVolume/tetrVol/CoeffQuality);
4519   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
4520   int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
4521   std::vector<int> aVec(SMDSEntity_Last);
4522   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
4523   if( IsQuadratic ) {
4524     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
4525     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
4526     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
4527   }
4528   else {
4529     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
4530     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
4531     aVec[SMDSEntity_Pyramid] = nbqua;
4532   }
4533   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
4534   aResMap.insert(std::make_pair(sm,aVec));
4535
4536   return true;
4537 }
4538
4539 bool GHS3DPlugin_GHS3D::importGMFMesh(const char* theGMFFileName, SMESH_Mesh& theMesh)
4540 {
4541   SMESH_MesherHelper* helper  = new SMESH_MesherHelper(theMesh );
4542 //   TopoDS_Shape theShape = theMesh.GetShapeToMesh();
4543   std::vector <const SMDS_MeshNode*> dummyNodeVector;
4544   std::map<const SMDS_MeshNode*,int> dummyNodeMap;
4545   std::map<std::vector<double>, std::string> dummyEnfVertGroup;
4546   std::vector<std::string> dummyElemGroup;
4547   std::set<std::string> dummyGroupsToRemove;
4548
4549   bool ok = readGMFFile(theGMFFileName,
4550                         this,
4551                         helper, theMesh.GetShapeToMesh(), dummyNodeVector, dummyNodeMap, dummyElemGroup, dummyElemGroup, dummyElemGroup, dummyGroupsToRemove);
4552   theMesh.GetMeshDS()->Modified();
4553   return ok;
4554 }
4555
4556 namespace
4557 {
4558   //================================================================================
4559   /*!
4560    * \brief Sub-mesh event listener setting enforced elements as soon as an enforced
4561    *        mesh is loaded
4562    */
4563   struct _EnforcedMeshRestorer : public SMESH_subMeshEventListener
4564   {
4565     _EnforcedMeshRestorer():
4566       SMESH_subMeshEventListener( /*isDeletable = */true, Name() )
4567     {}
4568
4569     //================================================================================
4570     /*!
4571      * \brief Returns an ID of listener
4572      */
4573     static const char* Name() { return "GHS3DPlugin_GHS3D::_EnforcedMeshRestorer"; }
4574
4575     //================================================================================
4576     /*!
4577      * \brief Treat events of the subMesh
4578      */
4579     void ProcessEvent(const int                       event,
4580                       const int                       eventType,
4581                       SMESH_subMesh*                  subMesh,
4582                       SMESH_subMeshEventListenerData* data,
4583                       const SMESH_Hypothesis*         hyp)
4584     {
4585       if ( SMESH_subMesh::SUBMESH_LOADED == event &&
4586            SMESH_subMesh::COMPUTE_EVENT  == eventType &&
4587            data &&
4588            !data->mySubMeshes.empty() )
4589       {
4590         // An enforced mesh (subMesh->_father) has been loaded from hdf file
4591         if ( GHS3DPlugin_Hypothesis* hyp = GetGHSHypothesis( data->mySubMeshes.front() ))
4592           hyp->RestoreEnfElemsByMeshes();
4593       }
4594     }
4595     //================================================================================
4596     /*!
4597      * \brief Returns GHS3DPlugin_Hypothesis used to compute a subMesh
4598      */
4599     static GHS3DPlugin_Hypothesis* GetGHSHypothesis( SMESH_subMesh* subMesh )
4600     {
4601       SMESH_HypoFilter ghsHypFilter( SMESH_HypoFilter::HasName( "GHS3D_Parameters" ));
4602       return (GHS3DPlugin_Hypothesis* )
4603         subMesh->GetFather()->GetHypothesis( subMesh->GetSubShape(),
4604                                              ghsHypFilter,
4605                                              /*visitAncestors=*/true);
4606     }
4607   };
4608
4609   //================================================================================
4610   /*!
4611    * \brief Sub-mesh event listener removing empty groups created due to "To make
4612    *        groups of domains".
4613    */
4614   struct _GroupsOfDomainsRemover : public SMESH_subMeshEventListener
4615   {
4616     _GroupsOfDomainsRemover():
4617       SMESH_subMeshEventListener( /*isDeletable = */true,
4618                                   "GHS3DPlugin_GHS3D::_GroupsOfDomainsRemover" ) {}
4619     /*!
4620      * \brief Treat events of the subMesh
4621      */
4622     void ProcessEvent(const int                       event,
4623                       const int                       eventType,
4624                       SMESH_subMesh*                  subMesh,
4625                       SMESH_subMeshEventListenerData* data,
4626                       const SMESH_Hypothesis*         hyp)
4627     {
4628       if (SMESH_subMesh::ALGO_EVENT == eventType &&
4629           !subMesh->GetAlgo() )
4630       {
4631         removeEmptyGroupsOfDomains( subMesh->GetFather(), /*notEmptyAsWell=*/true );
4632       }
4633     }
4634   };
4635 }
4636
4637 //================================================================================
4638 /*!
4639  * \brief Set an event listener to set enforced elements as soon as an enforced
4640  *        mesh is loaded
4641  */
4642 //================================================================================
4643
4644 void GHS3DPlugin_GHS3D::SubmeshRestored(SMESH_subMesh* subMesh)
4645 {
4646   if ( GHS3DPlugin_Hypothesis* hyp = _EnforcedMeshRestorer::GetGHSHypothesis( subMesh ))
4647   {
4648     GHS3DPlugin_Hypothesis::TGHS3DEnforcedMeshList enfMeshes = hyp->_GetEnforcedMeshes();
4649     GHS3DPlugin_Hypothesis::TGHS3DEnforcedMeshList::iterator it = enfMeshes.begin();
4650     for(;it != enfMeshes.end();++it) {
4651       GHS3DPlugin_Hypothesis::TGHS3DEnforcedMesh* enfMesh = *it;
4652       if ( SMESH_Mesh* mesh = GetMeshByPersistentID( enfMesh->persistID ))
4653       {
4654         SMESH_subMesh* smToListen = mesh->GetSubMesh( mesh->GetShapeToMesh() );
4655         // a listener set to smToListen will care of hypothesis stored in SMESH_EventListenerData
4656         subMesh->SetEventListener( new _EnforcedMeshRestorer(),
4657                                    SMESH_subMeshEventListenerData::MakeData( subMesh ),
4658                                    smToListen);
4659       }
4660     }
4661   }
4662 }
4663
4664 //================================================================================
4665 /*!
4666  * \brief Sets an event listener removing empty groups created due to "To make
4667  *        groups of domains".
4668  * \param subMesh - submesh where algo is set
4669  *
4670  * This method is called when a submesh gets HYP_OK algo_state.
4671  * After being set, event listener is notified on each event of a submesh.
4672  */
4673 //================================================================================
4674
4675 void GHS3DPlugin_GHS3D::SetEventListener(SMESH_subMesh* subMesh)
4676 {
4677   subMesh->SetEventListener( new _GroupsOfDomainsRemover(), 0, subMesh );
4678 }