src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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