src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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