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