src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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