src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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