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