src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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