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 #include <SMESH_subMeshEventListener.hxx>
  41 #include <SMESH_HypoFilter.hxx>
  42
  43 #include "SMDS_MeshElement.hxx"
  44 #include "SMDS_MeshNode.hxx"
  45 #include "SMDS_FaceOfNodes.hxx"
  46 #include "SMDS_VolumeOfNodes.hxx"
  47
  48 #include "SMESHDS_Group.hxx"
  49
  50 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  51 #include <StdMeshers_ViscousLayers.hxx>
  52
  53 #include <BRepAdaptor_Surface.hxx>
  54 #include <BRepBndLib.hxx>
  55 #include <BRepBuilderAPI_MakeVertex.hxx>
  56 #include <BRepClass3d_SolidClassifier.hxx>
  57 #include <BRepExtrema_DistShapeShape.hxx>
  58 #include <BRepGProp.hxx>
  59 #include <BRepTools.hxx>
  60 #include <BRep_Tool.hxx>
  61 #include <Bnd_Box.hxx>
  62 #include <GProp_GProps.hxx>
  63 #include <GeomAPI_ProjectPointOnSurf.hxx>
  64 #include <OSD_File.hxx>
  65 #include <Precision.hxx>
  66 #include <Quantity_Parameter.hxx>
  67 #include <Standard_ErrorHandler.hxx>
  68 #include <Standard_Failure.hxx>
  69 #include <Standard_ProgramError.hxx>
  70 #include <TopExp.hxx>
  71 #include <TopExp_Explorer.hxx>
  72 #include <TopTools_IndexedMapOfShape.hxx>
  73 #include <TopTools_ListIteratorOfListOfShape.hxx>
  74 #include <TopTools_MapOfShape.hxx>
  75 #include <TopoDS.hxx>
  76 #include <TopoDS_Shape.hxx>
  77 #include <TopoDS_Solid.hxx>
  78
  79 #include "utilities.h"
  80
  81 #ifdef WIN32
  82 #include <io.h>
  83 #else
  84 #include <sys/sysinfo.h>
  85 #endif
  86 #include <algorithm>
  87
  88 //#include <Standard_Stream.hxx>
  89
  90
  91 #define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
  92
  93 #ifdef _DEBUG_
  94 #define DUMP(txt) \
  95 //  std::cout << txt
  96 #else
  97 #define DUMP(txt)
  98 #endif
  99
 100 extern "C"
 101 {
 102 #ifndef WNT
 103 #include <unistd.h>
 104 #include <sys/mman.h>
 105 #endif
 106 #include <sys/stat.h>
 107 #include <fcntl.h>
 108 }
 109
 110 #define HOLE_ID -1
 111
 112 typedef const list<const SMDS_MeshFace*> TTriaList;
 113
 114 static void removeFile( const TCollection_AsciiString& fileName )
 115 {
 116   try {
 117     OSD_File( fileName ).Remove();
 118   }
 119   catch ( Standard_ProgramError ) {
 120     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
 121   }
 122 }
 123
 124 //=============================================================================
 125 /*!
 126  *
 127  */
 128 //=============================================================================
 129
 130 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
 131   : SMESH_3D_Algo(hypId, studyId, gen)
 132 {
 133   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
 134   _name = "GHS3D_3D";
 135   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
 136   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
 137   _iShape=0;
 138   _nbShape=0;
 139   _compatibleHypothesis.push_back( GHS3DPlugin_Hypothesis::GetHypType());
 140   _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() );
 141   _requireShape = false; // can work without shape_studyId
 142
 143   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 144   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 145   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 146
 147   MESSAGE("studyid = " << _studyId);
 148
 149   myStudy = NULL;
 150   myStudy = aStudyMgr->GetStudyByID(_studyId);
 151   if (myStudy)
 152     MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
 153
 154 #ifdef WITH_SMESH_CANCEL_COMPUTE
 155   _compute_canceled = false;
 156 #endif
 157 }
 158
 159 //=============================================================================
 160 /*!
 161  *
 162  */
 163 //=============================================================================
 164
 165 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
 166 {
 167   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
 168 }
 169
 170 //=============================================================================
 171 /*!
 172  *
 173  */
 174 //=============================================================================
 175
 176 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&         aMesh,
 177                                           const TopoDS_Shape& aShape,
 178                                           Hypothesis_Status&  aStatus )
 179 {
 180   aStatus = SMESH_Hypothesis::HYP_OK;
 181
 182   _hyp = 0;
 183   _viscousLayersHyp = 0;
 184   _keepFiles = false;
 185
 186   const list <const SMESHDS_Hypothesis * >& hyps =
 187     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 188   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 189   for ( ; h != hyps.end(); ++h )
 190   {
 191     if ( !_hyp )
 192       _hyp = dynamic_cast< const GHS3DPlugin_Hypothesis*> ( *h );
 193     if ( !_viscousLayersHyp )
 194       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 195   }
 196   if ( _hyp )
 197     _keepFiles = _hyp->GetKeepFiles();
 198
 199   return true;
 200 }
 201
 202
 203 //=======================================================================
 204 //function : entryToShape
 205 //purpose  :
 206 //=======================================================================
 207
 208 TopoDS_Shape GHS3DPlugin_GHS3D::entryToShape(std::string entry)
 209 {
 210   MESSAGE("GHS3DPlugin_GHS3D::entryToShape "<<entry );
 211   GEOM::GEOM_Object_var aGeomObj;
 212   TopoDS_Shape S = TopoDS_Shape();
 213   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 214   SALOMEDS::GenericAttribute_var anAttr;
 215
 216   if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
 217     SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 218     CORBA::String_var aVal = anIOR->Value();
 219     CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
 220     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 221   }
 222   if ( !aGeomObj->_is_nil() )
 223     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 224   return S;
 225 }
 226
 227 //=======================================================================
 228 //function : findShape
 229 //purpose  :
 230 //=======================================================================
 231
 232 static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[],
 233                               TopoDS_Shape        aShape,
 234                               const TopoDS_Shape  shape[],
 235                               double**            box,
 236                               const int           nShape,
 237                               TopAbs_State *      state = 0)
 238 {
 239   gp_XYZ aPnt(0,0,0);
 240   int j, iShape, nbNode = 4;
 241
 242   for ( j=0; j<nbNode; j++ ) {
 243     gp_XYZ p ( aNode[j]->X(), aNode[j]->Y(), aNode[j]->Z() );
 244     if ( aNode[j]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) {
 245       aPnt = p;
 246       break;
 247     }
 248     aPnt += p / nbNode;
 249   }
 250
 251   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 252   if (state) *state = SC.State();
 253   if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) {
 254     for (iShape = 0; iShape < nShape; iShape++) {
 255       aShape = shape[iShape];
 256       if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() ||
 257               aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() ||
 258               aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) {
 259         BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 260         if (state) *state = SC.State();
 261         if (SC.State() == TopAbs_IN)
 262           break;
 263       }
 264     }
 265   }
 266   return aShape;
 267 }
 268
 269 //=======================================================================
 270 //function : readMapIntLine
 271 //purpose  :
 272 //=======================================================================
 273
 274 static char* readMapIntLine(char* ptr, int tab[]) {
 275   long int intVal;
 276   std::cout << std::endl;
 277
 278   for ( int i=0; i<17; i++ ) {
 279     intVal = strtol(ptr, &ptr, 10);
 280     if ( i < 3 )
 281       tab[i] = intVal;
 282   }
 283   return ptr;
 284 }
 285
 286 //================================================================================
 287 /*!
 288  * \brief returns true if a triangle defined by the nodes is a temporary face on a
 289  * side facet of pyramid and defines sub-domian inside the pyramid
 290  */
 291 //================================================================================
 292
 293 static bool isTmpFace(const SMDS_MeshNode* node1,
 294                       const SMDS_MeshNode* node2,
 295                       const SMDS_MeshNode* node3)
 296 {
 297   // find a pyramid sharing the 3 nodes
 298   //const SMDS_MeshElement* pyram = 0;
 299   SMDS_ElemIteratorPtr vIt1 = node1->GetInverseElementIterator(SMDSAbs_Volume);
 300   while ( vIt1->more() )
 301   {
 302     const SMDS_MeshElement* pyram = vIt1->next();
 303     if ( pyram->NbCornerNodes() != 5 ) continue;
 304     int i2, i3;
 305     if ( (i2 = pyram->GetNodeIndex( node2 )) >= 0 &&
 306          (i3 = pyram->GetNodeIndex( node3 )) >= 0 )
 307     {
 308       // Triangle defines sub-domian inside the pyramid if it's
 309       // normal points out of the pyram
 310
 311       // make i2 and i3 hold indices of base nodes of the pyram while
 312       // keeping the nodes order in the triangle
 313       const int iApex = 4;
 314       if ( i2 == iApex )
 315         i2 = i3, i3 = pyram->GetNodeIndex( node1 );
 316       else if ( i3 == iApex )
 317         i3 = i2, i2 = pyram->GetNodeIndex( node1 );
 318
 319       int i3base = (i2+1) % 4; // next index after i2 within the pyramid base
 320       return ( i3base != i3 );
 321     }
 322   }
 323   return false;
 324 }
 325
 326 //=======================================================================
 327 //function : findShapeID
 328 //purpose  : find the solid corresponding to GHS3D sub-domain following
 329 //           the technique proposed in GHS3D manual (available within
 330 //           ghs3d installation) in chapter "B.4 Subdomain (sub-region) assignment".
 331 //           In brief: normal of the triangle defined by the given nodes
 332 //           points out of the domain it is associated to
 333 //=======================================================================
 334
 335 static int findShapeID(SMESH_Mesh&          mesh,
 336                        const SMDS_MeshNode* node1,
 337                        const SMDS_MeshNode* node2,
 338                        const SMDS_MeshNode* node3,
 339                        const bool           toMeshHoles)
 340 {
 341   const int invalidID = 0;
 342   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 343
 344   // face the nodes belong to
 345   const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
 346   if ( !face )
 347     return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
 348 #ifdef _DEBUG_
 349   std::cout << "bnd face " << face->GetID() << " - ";
 350 #endif
 351   // geom face the face assigned to
 352   SMESH_MeshEditor editor(&mesh);
 353   int geomFaceID = editor.FindShape( face );
 354   if ( !geomFaceID )
 355     return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
 356   TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
 357   if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
 358     return invalidID;
 359   TopoDS_Face geomFace = TopoDS::Face( shape );
 360
 361   // solids bounded by geom face
 362   TopTools_IndexedMapOfShape solids, shells;
 363   TopTools_ListIteratorOfListOfShape ansIt = mesh.GetAncestors(geomFace);
 364   for ( ; ansIt.More(); ansIt.Next() ) {
 365     switch ( ansIt.Value().ShapeType() ) {
 366     case TopAbs_SOLID:
 367       solids.Add( ansIt.Value() ); break;
 368     case TopAbs_SHELL:
 369       shells.Add( ansIt.Value() ); break;
 370     default:;
 371     }
 372   }
 373   // analyse found solids
 374   if ( solids.Extent() == 0 || shells.Extent() == 0)
 375     return invalidID;
 376
 377   const TopoDS_Solid& solid1 = TopoDS::Solid( solids(1) );
 378   if ( solids.Extent() == 1 )
 379   {
 380     if ( toMeshHoles )
 381       return meshDS->ShapeToIndex( solid1 );
 382
 383     //////////// UNCOMMENT AS SOON AS
 384     //////////// http://tracker.dev.opencascade.org/view.php?id=23129
 385     //////////// IS SOLVED
 386     // - Are we at a hole boundary face?
 387     // if ( shells(1).IsSame( BRepTools::OuterShell( solid1 )) )
 388     // { // - No, but maybe a hole is bound by two shapes? Does shells(1) touches another shell?
 389     //   bool touch = false;
 390     //   TopExp_Explorer eExp( shells(1), TopAbs_EDGE );
 391     //   // check if any edge of shells(1) belongs to another shell
 392     //   for ( ; eExp.More() && !touch; eExp.Next() ) {
 393     //     ansIt = mesh.GetAncestors( eExp.Current() );
 394     //     for ( ; ansIt.More() && !touch; ansIt.Next() ) {
 395     //       if ( ansIt.Value().ShapeType() == TopAbs_SHELL )
 396     //         touch = ( !ansIt.Value().IsSame( shells(1) ));
 397     //     }
 398     //   }
 399     //   if (!touch)
 400     //     return meshDS->ShapeToIndex( solid1 );
 401     // }
 402   }
 403   // find orientation of geom face within the first solid
 404   TopExp_Explorer fExp( solid1, TopAbs_FACE );
 405   for ( ; fExp.More(); fExp.Next() )
 406     if ( geomFace.IsSame( fExp.Current() )) {
 407       geomFace = TopoDS::Face( fExp.Current() );
 408       break;
 409     }
 410   if ( !fExp.More() )
 411     return invalidID; // face not found
 412
 413   // normale to triangle
 414   gp_Pnt node1Pnt ( node1->X(), node1->Y(), node1->Z() );
 415   gp_Pnt node2Pnt ( node2->X(), node2->Y(), node2->Z() );
 416   gp_Pnt node3Pnt ( node3->X(), node3->Y(), node3->Z() );
 417   gp_Vec vec12( node1Pnt, node2Pnt );
 418   gp_Vec vec13( node1Pnt, node3Pnt );
 419   gp_Vec meshNormal = vec12 ^ vec13;
 420   if ( meshNormal.SquareMagnitude() < DBL_MIN )
 421     return invalidID;
 422
 423   // get normale to geomFace at any node
 424   bool geomNormalOK = false;
 425   gp_Vec geomNormal;
 426   const SMDS_MeshNode* nodes[3] = { node1, node2, node3 };
 427   SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace );
 428   for ( int i = 0; !geomNormalOK && i < 3; ++i )
 429   {
 430     // find UV of i-th node on geomFace
 431     const SMDS_MeshNode* nNotOnSeamEdge = 0;
 432     if ( helper.IsSeamShape( nodes[i]->getshapeId() )) {
 433       if ( helper.IsSeamShape( nodes[(i+1)%3]->getshapeId() ))
 434         nNotOnSeamEdge = nodes[(i+2)%3];
 435       else
 436         nNotOnSeamEdge = nodes[(i+1)%3];
 437     }
 438     bool uvOK;
 439     gp_XY uv = helper.GetNodeUV( geomFace, nodes[i], nNotOnSeamEdge, &uvOK );
 440     // check that uv is correct
 441     if (uvOK) {
 442       double tol = 1e-6;
 443       TopoDS_Shape nodeShape = helper.GetSubShapeByNode( nodes[i], meshDS );
 444       if ( !nodeShape.IsNull() )
 445         switch ( nodeShape.ShapeType() )
 446         {
 447         case TopAbs_FACE:   tol = BRep_Tool::Tolerance( TopoDS::Face( nodeShape )); break;
 448         case TopAbs_EDGE:   tol = BRep_Tool::Tolerance( TopoDS::Edge( nodeShape )); break;
 449         case TopAbs_VERTEX: tol = BRep_Tool::Tolerance( TopoDS::Vertex( nodeShape )); break;
 450         default:;
 451         }
 452       gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
 453       BRepAdaptor_Surface surface( geomFace );
 454       uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
 455       if ( uvOK ) {
 456         // normale to geomFace at UV
 457         gp_Vec du, dv;
 458         surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
 459         geomNormal = du ^ dv;
 460         if ( geomFace.Orientation() == TopAbs_REVERSED )
 461           geomNormal.Reverse();
 462         geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
 463       }
 464     }
 465   }
 466   if ( !geomNormalOK)
 467     return invalidID;
 468
 469   // compare normals
 470   bool isReverse = ( meshNormal * geomNormal ) < 0;
 471   if ( !isReverse )
 472     return meshDS->ShapeToIndex( solid1 );
 473
 474   if ( solids.Extent() == 1 )
 475     return HOLE_ID; // we are inside a hole
 476   else
 477     return meshDS->ShapeToIndex( solids(2) );
 478 }
 479
 480 // //=======================================================================
 481 // //function : countShape
 482 // //purpose  :
 483 // //=======================================================================
 484 //
 485 // template < class Mesh, class Shape >
 486 // static int countShape( Mesh* mesh, Shape shape ) {
 487 //   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
 488 //   TopTools_MapOfShape mapShape;
 489 //   int nbShape = 0;
 490 //   for ( ; expShape.More(); expShape.Next() ) {
 491 //     if (mapShape.Add(expShape.Current())) {
 492 //       nbShape++;
 493 //     }
 494 //   }
 495 //   return nbShape;
 496 // }
 497 //
 498 // //=======================================================================
 499 // //function : getShape
 500 // //purpose  :
 501 // //=======================================================================
 502 //
 503 // template < class Mesh, class Shape, class Tab >
 504 // void getShape(Mesh* mesh, Shape shape, Tab *t_Shape) {
 505 //   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
 506 //   TopTools_MapOfShape mapShape;
 507 //   for ( int i=0; expShape.More(); expShape.Next() ) {
 508 //     if (mapShape.Add(expShape.Current())) {
 509 //       t_Shape[i] = expShape.Current();
 510 //       i++;
 511 //     }
 512 //   }
 513 //   return;
 514 // }
 515 //
 516 // // //=======================================================================
 517 // // //function : findEdgeID
 518 // // //purpose  :
 519 // // //=======================================================================
 520 //
 521 // static int findEdgeID(const SMDS_MeshNode* aNode,
 522 //                       const SMESHDS_Mesh*  theMesh,
 523 //                       const int            nEdge,
 524 //                       const TopoDS_Shape*  t_Edge) {
 525 //
 526 //   TopoDS_Shape aPntShape, foundEdge;
 527 //   TopoDS_Vertex aVertex;
 528 //   gp_Pnt aPnt( aNode->X(), aNode->Y(), aNode->Z() );
 529 //
 530 //   int foundInd, ind;
 531 //   double nearest = RealLast(), *t_Dist;
 532 //   double epsilon = Precision::Confusion();
 533 //
 534 //   t_Dist = new double[ nEdge ];
 535 //   aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape();
 536 //   aVertex   = TopoDS::Vertex( aPntShape );
 537 //
 538 //   for ( ind=0; ind < nEdge; ind++ ) {
 539 //     BRepExtrema_DistShapeShape aDistance ( aVertex, t_Edge[ind] );
 540 //     t_Dist[ind] = aDistance.Value();
 541 //     if ( t_Dist[ind] < nearest ) {
 542 //       nearest   = t_Dist[ind];
 543 //       foundEdge = t_Edge[ind];
 544 //       foundInd  = ind;
 545 //       if ( nearest < epsilon )
 546 //         ind = nEdge;
 547 //     }
 548 //   }
 549 //
 550 //   delete [] t_Dist;
 551 //   return theMesh->ShapeToIndex( foundEdge );
 552 // }
 553 //
 554 //
 555 // // =======================================================================
 556 // // function : readGMFFile
 557 // // purpose  : read GMF file with geometry associated to mesh
 558 // // =======================================================================
 559 //
 560 // static bool readGMFFile(const int                       fileOpen,
 561 //                         const char*                     theFileName,
 562 //                         SMESH_Mesh&                     theMesh,
 563 //                         const int                       nbShape,
 564 //                         const TopoDS_Shape*             tabShape,
 565 //                         double**                        tabBox,
 566 //                         map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
 567 //                         bool                            toMeshHoles,
 568 //                         int                             nbEnforcedVertices,
 569 //                         int                             nbEnforcedNodes)
 570 // {
 571 //   TopoDS_Shape aShape;
 572 //   TopoDS_Vertex aVertex;
 573 //   SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
 574 //   int nbElem = 0, nbRef = 0, IdShapeRef = 1;
 575 //   int *tabID;
 576 //   int aGMFNodeID = 0;
 577 //   int compoundID =
 578 //     nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
 579 //   int tetraShapeID = compoundID;
 580 //   double epsilon = Precision::Confusion();
 581 //   int *nodeAssigne, *GMFNodeAssigne;
 582 //   SMDS_MeshNode** GMFNode;
 583 //   TopoDS_Shape *tabCorner, *tabEdge;
 584 //   std::map <GmfKwdCod,int> tabRef;
 585 //
 586 //
 587 //   int ver, dim;
 588 //   MESSAGE("Read " << theFileName << " file");
 589 //   int InpMsh = GmfOpenMesh(theFileName, GmfRead, &ver, &dim);
 590 //   if (!InpMsh)
 591 //     return false;
 592 //
 593 //   // ===========================
 594 //   // Fill the tabID array: BEGIN
 595 //   // ===========================
 596 //
 597 //   /*
 598 //   The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
 599 //   */
 600 //   Kernel_Utils::Localizer loc;
 601 //   struct stat status;
 602 //   size_t      length;
 603 //
 604 //   char *ptr, *mapPtr;
 605 //   char *tetraPtr;
 606 //   int *tab = new int[3];
 607 //
 608 //   // Read the file state
 609 //   fstat(fileOpen, &status);
 610 //   length   = status.st_size;
 611 //
 612 //   // Mapping the result file into memory
 613 // #ifdef WNT
 614 //   HANDLE fd = CreateFile(theFileName, GENERIC_READ, FILE_SHARE_READ,
 615 //                          NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
 616 //   HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
 617 //                                         0, (DWORD)length, NULL);
 618 //   ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
 619 // #else
 620 //   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
 621 // #endif
 622 //   mapPtr = ptr;
 623 //
 624 //   ptr      = readMapIntLine(ptr, tab);
 625 //   tetraPtr = ptr;
 626 //
 627 //   nbElem            = tab[0];
 628 //   int nbNodes       = tab[1];
 629 //
 630 //   for (int i=0; i < 4*nbElem; i++)
 631 //     strtol(ptr, &ptr, 10);
 632 //
 633 //   for (int iNode=1; iNode <= nbNodes; iNode++)
 634 //     for (int iCoor=0; iCoor < 3; iCoor++)
 635 //       strtod(ptr, &ptr);
 636 //
 637 //
 638 //   // Reading the number of triangles which corresponds to the number of sub-domains
 639 //   int nbTriangle = strtol(ptr, &ptr, 10);
 640 //
 641 //
 642 //   // The keyword does not exist yet => to update when it is created
 643 // //   int nbSubdomains = GmfStatKwd(InpMsh, GmfSubdomain);
 644 // //   int id_tri[3];
 645 //
 646 //
 647 //   tabID = new int[nbTriangle];
 648 //   for (int i=0; i < nbTriangle; i++) {
 649 //     tabID[i] = 0;
 650 //     int nodeId1, nodeId2, nodeId3;
 651 //     // find the solid corresponding to GHS3D sub-domain following
 652 //     // the technique proposed in GHS3D manual in chapter
 653 //     // "B.4 Subdomain (sub-region) assignment"
 654 //
 655 //     nodeId1 = strtol(ptr, &ptr, 10);
 656 //     nodeId2 = strtol(ptr, &ptr, 10);
 657 //     nodeId3 = strtol(ptr, &ptr, 10);
 658 //
 659 // //   // The keyword does not exist yet => to update when it is created
 660 // //     GmfGetLin(InpMsh, GmfSubdomain, &id_tri[0], &id_tri[1], &id_tri[2]);
 661 // //     nodeId1 = id_tri[0];
 662 // //     nodeId2 = id_tri[1];
 663 // //     nodeId3 = id_tri[2];
 664 //
 665 //     if ( nbTriangle > 1 ) {
 666 //       // get the nodes indices
 667 //       const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
 668 //       const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
 669 //       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
 670 //       try {
 671 //         OCC_CATCH_SIGNALS;
 672 //         tabID[i] = findShapeID( theMesh, n1, n2, n3, toMeshHoles );
 673 //         // -- 0020330: Pb with ghs3d as a submesh
 674 //         // check that found shape is to be meshed
 675 //         if ( tabID[i] > 0 ) {
 676 //           const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] );
 677 //           bool isToBeMeshed = false;
 678 //           for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS )
 679 //             isToBeMeshed = foundShape.IsSame( tabShape[ iS ]);
 680 //           if ( !isToBeMeshed )
 681 //             tabID[i] = HOLE_ID;
 682 //         }
 683 //         // END -- 0020330: Pb with ghs3d as a submesh
 684 // #ifdef _DEBUG_
 685 //         std::cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << std::endl;
 686 // #endif
 687 //       }
 688 //       catch ( Standard_Failure & ex)
 689 //       {
 690 // #ifdef _DEBUG_
 691 //         std::cout << i+1 << " subdomain: Exception caugt: " << ex.GetMessageString() << std::endl;
 692 // #endif
 693 //       }
 694 //       catch (...) {
 695 // #ifdef _DEBUG_
 696 //         std::cout << i+1 << " subdomain: unknown exception caught " << std::endl;
 697 // #endif
 698 //       }
 699 //     }
 700 //   }
 701 //
 702 //   // ===========================
 703 //   // Fill the tabID array: END
 704 //   // ===========================
 705 //
 706 //
 707 //   tabRef[GmfVertices]       = 3;
 708 //   tabRef[GmfCorners]        = 1;
 709 //   tabRef[GmfEdges]          = 2;
 710 //   tabRef[GmfRidges]         = 1;
 711 //   tabRef[GmfTriangles]      = 3;
 712 // //   tabRef[GmfQuadrilaterals] = 4;
 713 //   tabRef[GmfTetrahedra]     = 4;
 714 // //   tabRef[GmfHexahedra]      = 8;
 715 //
 716 //   SMDS_NodeIteratorPtr itOnGMFInputNode = theMeshDS->nodesIterator();
 717 //   while ( itOnGMFInputNode->more() )
 718 //     theMeshDS->RemoveNode( itOnGMFInputNode->next() );
 719 //
 720 //
 721 //   int nbVertices = GmfStatKwd(InpMsh, GmfVertices);
 722 //   int nbCorners = max(countShape( theMeshDS, TopAbs_VERTEX ) , GmfStatKwd(InpMsh, GmfCorners));
 723 //   int nbShapeEdge = countShape( theMeshDS, TopAbs_EDGE );
 724 //
 725 //   tabCorner       = new TopoDS_Shape[ nbCorners ];
 726 //   tabEdge         = new TopoDS_Shape[ nbShapeEdge ];
 727 //   nodeAssigne     = new int[ nbVertices + 1 ];
 728 //   GMFNodeAssigne  = new int[ nbVertices + 1 ];
 729 //   GMFNode         = new SMDS_MeshNode*[ nbVertices + 1 ];
 730 //
 731 //   getShape(theMeshDS, TopAbs_VERTEX, tabCorner);
 732 //   getShape(theMeshDS, TopAbs_EDGE,   tabEdge);
 733 //
 734 //   std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
 735 //   for ( ; it != tabRef.end() ; ++it)
 736 //   {
 737 // //     int dummy;
 738 //     GmfKwdCod token = it->first;
 739 //     nbRef    = it->second;
 740 //
 741 //     nbElem = GmfStatKwd(InpMsh, token);
 742 //     if (nbElem > 0) {
 743 //       GmfGotoKwd(InpMsh, token);
 744 //       std::cout << "Read " << nbElem;
 745 //     }
 746 //     else
 747 //       continue;
 748 //
 749 //     int id[nbElem*tabRef[token]];
 750 //     int ghs3dShapeID[nbElem];
 751 //
 752 //     if (token == GmfVertices) {
 753 //       std::cout << " vertices" << std::endl;
 754 //       int aGMFID;
 755 //
 756 //       float VerTab_f[nbElem][3];
 757 //       double VerTab_d[nbElem][3];
 758 //       SMDS_MeshNode * aGMFNode;
 759 //
 760 //       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
 761 //         aGMFID = iElem + 1;
 762 //         if (ver == GmfFloat) {
 763 //           GmfGetLin(InpMsh, token, &VerTab_f[nbElem][0], &VerTab_f[nbElem][1], &VerTab_f[nbElem][2], &ghs3dShapeID[iElem]);
 764 //           aGMFNode = theMeshDS->AddNode(VerTab_f[nbElem][0], VerTab_f[nbElem][1], VerTab_f[nbElem][2]);
 765 //         }
 766 //         else {
 767 //           GmfGetLin(InpMsh, token, &VerTab_d[nbElem][0], &VerTab_d[nbElem][1], &VerTab_d[nbElem][2], &ghs3dShapeID[iElem]);
 768 //           aGMFNode = theMeshDS->AddNode(VerTab_d[nbElem][0], VerTab_d[nbElem][1], VerTab_d[nbElem][2]);
 769 //         }
 770 //         GMFNode[ aGMFID ] = aGMFNode;
 771 //         nodeAssigne[ aGMFID ] = 0;
 772 //         GMFNodeAssigne[ aGMFID ] = 0;
 773 //       }
 774 //     }
 775 //     else if (token == GmfCorners && nbElem > 0) {
 776 //       std::cout << " corners" << std::endl;
 777 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 778 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
 779 //     }
 780 //     else if (token == GmfRidges && nbElem > 0) {
 781 //       std::cout << " ridges" << std::endl;
 782 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 783 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
 784 //     }
 785 //     else if (token == GmfEdges && nbElem > 0) {
 786 //       std::cout << " edges" << std::endl;
 787 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 788 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &ghs3dShapeID[iElem]);
 789 //     }
 790 //     else if (token == GmfTriangles && nbElem > 0) {
 791 //       std::cout << " triangles" << std::endl;
 792 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 793 //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &ghs3dShapeID[iElem]);
 794 //     }
 795 // //     else if (token == GmfQuadrilaterals && nbElem > 0) {
 796 // //       std::cout << " Quadrilaterals" << std::endl;
 797 // //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 798 // //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &ghs3dShapeID[iElem]);
 799 // //     }
 800 //     else if (token == GmfTetrahedra && nbElem > 0) {
 801 //       std::cout << " Tetrahedra" << std::endl;
 802 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 803 //         GmfGetLin(InpMsh, token,
 804 //                   &id[iElem*tabRef[token]],
 805 //                   &id[iElem*tabRef[token]+1],
 806 //                   &id[iElem*tabRef[token]+2],
 807 //                   &id[iElem*tabRef[token]+3],
 808 //                   &ghs3dShapeID[iElem]);
 809 //     }
 810 // //     else if (token == GmfHexahedra && nbElem > 0) {
 811 // //       std::cout << " Hexahedra" << std::endl;
 812 // //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 813 // //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
 814 // //                   &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &ghs3dShapeID[iElem]);
 815 // //     }
 816 //
 817 //     switch (token) {
 818 //     case GmfCorners:
 819 //     case GmfRidges:
 820 //     case GmfEdges:
 821 //     case GmfTriangles:
 822 // //     case GmfQuadrilaterals:
 823 //     case GmfTetrahedra:
 824 // //     case GmfHexahedra:
 825 //     {
 826 //       int nodeDim, shapeID, *nodeID;
 827 //       const SMDS_MeshNode** node;
 828 // //       std::vector< SMDS_MeshNode* > enfNode( nbRef );
 829 //       SMDS_MeshElement * aGMFElement;
 830 //
 831 //       node    = new const SMDS_MeshNode*[nbRef];
 832 //       nodeID  = new int[ nbRef ];
 833 //
 834 //       for ( int iElem = 0; iElem < nbElem; iElem++ )
 835 //       {
 836 //         for ( int iRef = 0; iRef < nbRef; iRef++ )
 837 //         {
 838 //           aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
 839 //           node  [ iRef ] = GMFNode[ aGMFNodeID ];
 840 //           nodeID[ iRef ] = aGMFNodeID;
 841 //         }
 842 //
 843 //         switch (token)
 844 //         {
 845 //         case GmfCorners: {
 846 //           nodeDim = 1;
 847 //           gp_Pnt GMFPnt ( node[0]->X(), node[0]->Y(), node[0]->Z() );
 848 //           for ( int i=0; i<nbElem; i++ ) {
 849 //             aVertex = TopoDS::Vertex( tabCorner[i] );
 850 //             gp_Pnt aPnt = BRep_Tool::Pnt( aVertex );
 851 //             if ( aPnt.Distance( GMFPnt ) < epsilon )
 852 //               break;
 853 //           }
 854 //           break;
 855 //         }
 856 //         case GmfEdges: {
 857 //           nodeDim = 2;
 858 //           aGMFElement = theMeshDS->AddEdge( node[0], node[1] );
 859 //           int iNode = 1;
 860 //           if ( GMFNodeAssigne[ nodeID[0] ] == 0 || GMFNodeAssigne[ nodeID[0] ] == 2 )
 861 //             iNode = 0;
 862 //           shapeID = findEdgeID( node[iNode], theMeshDS, nbShapeEdge, tabEdge );
 863 //           break;
 864 //         }
 865 //         case GmfRidges:
 866 //           break;
 867 //         case GmfTriangles: {
 868 //           nodeDim = 3;
 869 //           aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2]);
 870 //           shapeID = -1;
 871 //           break;
 872 //         }
 873 // //         case GmfQuadrilaterals: {
 874 // //           nodeDim = 4;
 875 // //           aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2], node[3] );
 876 // //           shapeID = -1;
 877 // //           break;
 878 // //         }
 879 //         case GmfTetrahedra: {
 880 //
 881 //           // IN WORK
 882 //           TopoDS_Shape aSolid;
 883 //           // We always run GHS3D with "to mesh holes"==TRUE but we must not create
 884 //           // tetras within holes depending on hypo option,
 885 //           // so we first check if aTet is inside a hole and then create it
 886 //           if ( nbTriangle > 1 ) {
 887 //             tetraShapeID = HOLE_ID; // negative tetraShapeID means not to create tetras if !toMeshHoles
 888 //             int aGhs3dShapeID = ghs3dShapeID[iElem] - IdShapeRef;
 889 //             if ( tabID[ aGhs3dShapeID ] == 0 ) {
 890 //               TopAbs_State state;
 891 //               aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
 892 //               if ( toMeshHoles || state == TopAbs_IN )
 893 //                 tetraShapeID = theMeshDS->ShapeToIndex( aSolid );
 894 //               tabID[ aGhs3dShapeID ] = tetraShapeID;
 895 //             }
 896 //             else
 897 //               tetraShapeID = tabID[ aGhs3dShapeID ];
 898 //           }
 899 //           else if ( nbShape > 1 ) {
 900 //             // Case where nbTriangle == 1 while nbShape == 2 encountered
 901 //             // with compound of 2 boxes and "To mesh holes"==False,
 902 //             // so there are no subdomains specified for each tetrahedron.
 903 //             // Try to guess a solid by a node already bound to shape
 904 //             tetraShapeID = 0;
 905 //             for ( int i=0; i<4 && tetraShapeID==0; i++ ) {
 906 //               if ( nodeAssigne[ nodeID[i] ] == 1 &&
 907 //                   node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
 908 //                   node[i]->getshapeId() > 1 )
 909 //               {
 910 //                 tetraShapeID = node[i]->getshapeId();
 911 //               }
 912 //             }
 913 //             if ( tetraShapeID==0 ) {
 914 //               aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
 915 //               tetraShapeID = theMeshDS->ShapeToIndex( aSolid );
 916 //             }
 917 //           }
 918 //           // set new nodes and tetrahedron onto the shape
 919 //           for ( int i=0; i<4; i++ ) {
 920 //             if ( nodeAssigne[ nodeID[i] ] == 0 ) {
 921 //               if ( tetraShapeID != HOLE_ID )
 922 //                 theMeshDS->SetNodeInVolume( node[i], tetraShapeID );
 923 //               nodeAssigne[ nodeID[i] ] = tetraShapeID;
 924 //             }
 925 //           }
 926 //           if ( toMeshHoles || tetraShapeID != HOLE_ID ) {
 927 //             aGMFElement = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
 928 //             theMeshDS->SetMeshElementOnShape( aGMFElement, tetraShapeID );
 929 //           }
 930 //
 931 //           // IN WORK
 932 //
 933 //           nodeDim = 5;
 934 //           break;
 935 //         }
 936 // //         case GmfHexahedra: {
 937 // //           nodeDim = 6;
 938 // //           aGMFElement = theMeshDS->AddVolume( node[0], node[3], node[2], node[1],
 939 // //                                             node[4], node[7], node[6], node[5] );
 940 // //           break;
 941 // //         }
 942 //         default: continue;
 943 //         }
 944 //         if (token != GmfRidges)
 945 //         {
 946 //           for ( int i=0; i<nbRef; i++ ) {
 947 //               if ( GMFNodeAssigne[ nodeID[i] ] == 0 ) {
 948 //                 if      ( token == GmfCorners )   theMeshDS->SetNodeOnVertex( node[0], aVertex );
 949 //                 else if ( token == GmfEdges )     theMeshDS->SetNodeOnEdge( node[i], shapeID );
 950 //                 else if ( token == GmfTriangles ) theMeshDS->SetNodeOnFace( node[i], shapeID );
 951 //                 GMFNodeAssigne[ nodeID[i] ] = nodeDim;
 952 //               }
 953 //             }
 954 //             if ( token != "Corners" )
 955 //               theMeshDS->SetMeshElementOnShape( aGMFElement, shapeID );
 956 //         }
 957 //       } // for
 958 //
 959 //       if ( !toMeshHoles ) {
 960 //         map <int,const SMDS_MeshNode*>::iterator itOnNode = theGhs3dIdToNodeMap.find( nbVertices-(nbEnforcedVertices+nbEnforcedNodes) );
 961 //         for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
 962 //           if ( nodeAssigne[ itOnNode->first ] == HOLE_ID )
 963 //             theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
 964 //         }
 965 //       }
 966 //
 967 //       delete [] node;
 968 //       delete [] nodeID;
 969 //       break;
 970 //       } // case GmfTetrahedra
 971 //     } // switch(token)
 972 //   } // for
 973 //   cout << std::endl;
 974 //
 975 // #ifdef WNT
 976 //   UnmapViewOfFile(mapPtr);
 977 //   CloseHandle(hMapObject);
 978 //   CloseHandle(fd);
 979 // #else
 980 //   munmap(mapPtr, length);
 981 // #endif
 982 //   close(fileOpen);
 983 //
 984 //   delete [] tabID;
 985 //   delete [] tabCorner;
 986 //   delete [] tabEdge;
 987 //   delete [] nodeAssigne;
 988 //   delete [] GMFNodeAssigne;
 989 //   delete [] GMFNode;
 990 //
 991 //   return true;
 992 // }
 993
 994
 995 //=======================================================================
 996 //function : addElemInMeshGroup
 997 //purpose  : Update or create groups in mesh
 998 //=======================================================================
 999
1000 static void addElemInMeshGroup(SMESH_Mesh*             theMesh,
1001                                const SMDS_MeshElement* anElem,
1002                                std::string&            groupName,
1003                                std::set<std::string>&  groupsToRemove)
1004 {
1005   if ( !anElem ) return; // issue 0021776
1006
1007   bool groupDone = false;
1008   SMESH_Mesh::GroupIteratorPtr grIt = theMesh->GetGroups();
1009   while (grIt->more()) {
1010     SMESH_Group * group = grIt->next();
1011     if ( !group ) continue;
1012     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
1013     if ( !groupDS ) continue;
1014     if ( groupDS->GetType()==anElem->GetType() &&groupName.compare(group->GetName())==0) {
1015       SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
1016       aGroupDS->SMDSGroup().Add(anElem);
1017       groupDone = true;
1018 //       MESSAGE("Successfully added enforced element to existing group " << groupName);
1019       break;
1020     }
1021   }
1022
1023   if (!groupDone)
1024   {
1025     int groupId;
1026     SMESH_Group* aGroup = theMesh->AddGroup(anElem->GetType(), groupName.c_str(), groupId);
1027     aGroup->SetName( groupName.c_str() );
1028     SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
1029     aGroupDS->SMDSGroup().Add(anElem);
1030 //     MESSAGE("Successfully created enforced vertex group " << groupName);
1031     groupDone = true;
1032   }
1033   if (!groupDone)
1034     throw SALOME_Exception(LOCALIZED("A given element was not added to a group"));
1035 }
1036
1037
1038 //=======================================================================
1039 //function : updateMeshGroups
1040 //purpose  : Update or create groups in mesh
1041 //=======================================================================
1042
1043 static void updateMeshGroups(SMESH_Mesh* theMesh, std::set<std::string> groupsToRemove)
1044 {
1045   SMESH_Mesh::GroupIteratorPtr grIt = theMesh->GetGroups();
1046   while (grIt->more()) {
1047     SMESH_Group * group = grIt->next();
1048     if ( !group ) continue;
1049     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
1050     if ( !groupDS ) continue;
1051     std::string currentGroupName = (string)group->GetName();
1052     if (groupDS->IsEmpty() && groupsToRemove.find(currentGroupName) != groupsToRemove.end()) {
1053       // Previous group created by enforced elements
1054       MESSAGE("Delete previous group created by removed enforced elements: " << group->GetName())
1055       theMesh->RemoveGroup(groupDS->GetID());
1056     }
1057   }
1058 }
1059
1060 //=======================================================================
1061 //function : readGMFFile
1062 //purpose  : read GMF file w/o geometry associated to mesh
1063 //=======================================================================
1064
1065 static bool readGMFFile(const char*                     theFile,
1066 #ifdef WITH_SMESH_CANCEL_COMPUTE
1067                         GHS3DPlugin_GHS3D*              theAlgo,
1068 #endif
1069                         SMESH_MesherHelper*             theHelper,
1070                         TopoDS_Shape                    theSolid,
1071                         vector <const SMDS_MeshNode*> & theNodeByGhs3dId,
1072                         map<const SMDS_MeshNode*,int> & theNodeToGhs3dIdMap,
1073                         std::vector<std::string> &      aNodeGroupByGhs3dId,
1074                         std::vector<std::string> &      anEdgeGroupByGhs3dId,
1075                         std::vector<std::string> &      aFaceGroupByGhs3dId,
1076                         std::set<std::string> &         groupsToRemove
1077                        )
1078 {
1079   std::string tmpStr;
1080   SMESHDS_Mesh* theMeshDS = theHelper->GetMeshDS();
1081
1082   int nbInitialNodes = theNodeByGhs3dId.size();
1083   int nbMeshNodes = theMeshDS->NbNodes();
1084
1085   const bool isQuadMesh =
1086     theHelper->GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
1087     theHelper->GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
1088     theHelper->GetMesh()->NbVolumes( ORDER_QUADRATIC );
1089
1090 #ifdef _DEBUG_
1091   std::cout << "theNodeByGhs3dId.size(): " << nbInitialNodes << std::endl;
1092   std::cout << "theHelper->GetMesh()->NbNodes(): " << nbMeshNodes << std::endl;
1093   std::cout << "isQuadMesh: " << isQuadMesh << std::endl;
1094 #endif
1095
1096   if (theHelper->GetSubShapeID() != 0)
1097     theHelper->IsQuadraticSubMesh( theHelper->GetSubShape() );
1098
1099   // ---------------------------------
1100   // Read generated elements and nodes
1101   // ---------------------------------
1102
1103   int nbElem = 0, nbRef = 0;
1104   int aGMFNodeID = 0/*, shapeID*/;
1105   //int *nodeAssigne;
1106   const SMDS_MeshNode** GMFNode;
1107 #ifdef _DEBUG_
1108   std::map<int, std::set<int> > subdomainId2tetraId;
1109 #endif
1110   std::map <GmfKwdCod,int> tabRef;
1111
1112   tabRef[GmfVertices]       = 3; // for new nodes and enforced nodes
1113   tabRef[GmfCorners]        = 1;
1114   tabRef[GmfEdges]          = 2; // for enforced edges
1115   tabRef[GmfRidges]         = 1;
1116   tabRef[GmfTriangles]      = 3; // for enforced faces
1117   tabRef[GmfQuadrilaterals] = 4;
1118   tabRef[GmfTetrahedra]     = 4; // for new tetras
1119   tabRef[GmfHexahedra]      = 8;
1120
1121   int ver, dim;
1122   MESSAGE("Read " << theFile << " file");
1123   int InpMsh = GmfOpenMesh(theFile, GmfRead, &ver, &dim);
1124   if (!InpMsh)
1125     return false;
1126   MESSAGE("Done ");
1127
1128   // Issue 0020682. Avoid creating nodes and tetras at place where
1129   // volumic elements already exist
1130   SMESH_ElementSearcher* elemSearcher = 0;
1131   vector< const SMDS_MeshElement* > foundVolumes;
1132   if ( theHelper->GetMesh()->NbVolumes() > 0 )
1133     elemSearcher = SMESH_MeshEditor( theHelper->GetMesh() ).GetElementSearcher();
1134
1135   int nbVertices = GmfStatKwd(InpMsh, GmfVertices) - nbInitialNodes;
1136   GMFNode = new const SMDS_MeshNode*[ nbVertices + 1 ];
1137   //nodeAssigne = new int[ nbVertices + 1 ];
1138
1139   std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
1140   for ( ; it != tabRef.end() ; ++it)
1141   {
1142 #ifdef WITH_SMESH_CANCEL_COMPUTE
1143     if(theAlgo->computeCanceled()) {
1144       GmfCloseMesh(InpMsh);
1145       delete [] GMFNode;
1146       //delete [] nodeAssigne;
1147       return false;
1148     }
1149 #endif
1150     int dummy;
1151     GmfKwdCod token = it->first;
1152     nbRef    = it->second;
1153
1154     nbElem = GmfStatKwd(InpMsh, token);
1155     if (nbElem > 0) {
1156       GmfGotoKwd(InpMsh, token);
1157       std::cout << "Read " << nbElem;
1158     }
1159     else
1160       continue;
1161
1162     std::vector<int> id (nbElem*tabRef[token]); // node ids
1163
1164     if (token == GmfVertices) {
1165       (nbElem <= 1) ? tmpStr = " vertex" : tmpStr = " vertices";
1166 //       std::cout << nbInitialNodes << " from input mesh " << std::endl;
1167
1168       // Remove orphan nodes from previous enforced mesh which was cleared
1169 //       if ( nbElem < nbMeshNodes ) {
1170 //         const SMDS_MeshNode* node;
1171 //         SMDS_NodeIteratorPtr nodeIt = theMeshDS->nodesIterator();
1172 //         while ( nodeIt->more() )
1173 //         {
1174 //           node = nodeIt->next();
1175 //           if (theNodeToGhs3dIdMap.find(node) != theNodeToGhs3dIdMap.end())
1176 //             theMeshDS->RemoveNode(node);
1177 //         }
1178 //       }
1179
1180
1181       int aGMFID;
1182
1183       float VerTab_f[3];
1184       double x, y, z;
1185       const SMDS_MeshNode * aGMFNode;
1186
1187       //shapeID = theMeshDS->ShapeToIndex( theSolid );
1188       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
1189 #ifdef WITH_SMESH_CANCEL_COMPUTE
1190         if(theAlgo->computeCanceled()) {
1191           GmfCloseMesh(InpMsh);
1192           delete [] GMFNode;
1193           //delete [] nodeAssigne;
1194           return false;
1195         }
1196 #endif
1197         if (ver == GmfFloat) {
1198           GmfGetLin(InpMsh, token, &VerTab_f[0], &VerTab_f[1], &VerTab_f[2], &dummy);
1199           x = VerTab_f[0];
1200           y = VerTab_f[1];
1201           z = VerTab_f[2];
1202         }
1203         else {
1204           GmfGetLin(InpMsh, token, &x, &y, &z, &dummy);
1205         }
1206         if (iElem >= nbInitialNodes) {
1207           if ( elemSearcher &&
1208                 elemSearcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_Volume, foundVolumes))
1209             aGMFNode = 0;
1210           else
1211             aGMFNode = theHelper->AddNode(x, y, z);
1212
1213           aGMFID = iElem -nbInitialNodes +1;
1214           GMFNode[ aGMFID ] = aGMFNode;
1215           //nodeAssigne[ aGMFID ] = 0;
1216           if (aGMFID-1 < aNodeGroupByGhs3dId.size() && !aNodeGroupByGhs3dId.at(aGMFID-1).empty())
1217             addElemInMeshGroup(theHelper->GetMesh(), aGMFNode, aNodeGroupByGhs3dId.at(aGMFID-1), groupsToRemove);
1218         }
1219       }
1220     }
1221     else if (token == GmfCorners && nbElem > 0) {
1222       (nbElem <= 1) ? tmpStr = " corner" : tmpStr = " corners";
1223       for ( int iElem = 0; iElem < nbElem; iElem++ )
1224         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
1225     }
1226     else if (token == GmfRidges && nbElem > 0) {
1227       (nbElem <= 1) ? tmpStr = " ridge" : tmpStr = " ridges";
1228       for ( int iElem = 0; iElem < nbElem; iElem++ )
1229         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
1230     }
1231     else if (token == GmfEdges && nbElem > 0) {
1232       (nbElem <= 1) ? tmpStr = " edge" : tmpStr = " edges";
1233       for ( int iElem = 0; iElem < nbElem; iElem++ )
1234         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &dummy);
1235     }
1236     else if (token == GmfTriangles && nbElem > 0) {
1237       (nbElem <= 1) ? tmpStr = " triangle" : tmpStr = " triangles";
1238       for ( int iElem = 0; iElem < nbElem; iElem++ )
1239         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &dummy);
1240     }
1241     else if (token == GmfQuadrilaterals && nbElem > 0) {
1242       (nbElem <= 1) ? tmpStr = " Quadrilateral" : tmpStr = " Quadrilaterals";
1243       for ( int iElem = 0; iElem < nbElem; iElem++ )
1244         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &dummy);
1245     }
1246     else if (token == GmfTetrahedra && nbElem > 0) {
1247       (nbElem <= 1) ? tmpStr = " Tetrahedron" : tmpStr = " Tetrahedra";
1248       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
1249         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &dummy);
1250 #ifdef _DEBUG_
1251         subdomainId2tetraId[dummy].insert(iElem+1);
1252 //         MESSAGE("subdomainId2tetraId["<<dummy<<"].insert("<<iElem+1<<")");
1253 #endif
1254       }
1255     }
1256     else if (token == GmfHexahedra && nbElem > 0) {
1257       (nbElem <= 1) ? tmpStr = " Hexahedron" : tmpStr = " Hexahedra";
1258       for ( int iElem = 0; iElem < nbElem; iElem++ )
1259         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
1260                   &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &dummy);
1261     }
1262     std::cout << tmpStr << std::endl;
1263     std::cout << std::endl;
1264
1265     switch (token) {
1266     case GmfCorners:
1267     case GmfRidges:
1268     case GmfEdges:
1269     case GmfTriangles:
1270     case GmfQuadrilaterals:
1271     case GmfTetrahedra:
1272     case GmfHexahedra:
1273     {
1274       std::vector< const SMDS_MeshNode* > node( nbRef );
1275       std::vector< int >          nodeID( nbRef );
1276       std::vector< SMDS_MeshNode* > enfNode( nbRef );
1277       const SMDS_MeshElement* aCreatedElem;
1278
1279       for ( int iElem = 0; iElem < nbElem; iElem++ )
1280       {
1281 #ifdef WITH_SMESH_CANCEL_COMPUTE
1282         if(theAlgo->computeCanceled()) {
1283           GmfCloseMesh(InpMsh);
1284           delete [] GMFNode;
1285           //delete [] nodeAssigne;
1286           return false;
1287         }
1288 #endif
1289         // Check if elem is already in input mesh. If yes => skip
1290         bool fullyCreatedElement = false; // if at least one of the nodes was created
1291         for ( int iRef = 0; iRef < nbRef; iRef++ )
1292         {
1293           aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
1294           if (aGMFNodeID <= nbInitialNodes) // input nodes
1295           {
1296             aGMFNodeID--;
1297             node[ iRef ] = theNodeByGhs3dId[aGMFNodeID];
1298           }
1299           else
1300           {
1301             fullyCreatedElement = true;
1302             aGMFNodeID -= nbInitialNodes;
1303             nodeID[ iRef ] = aGMFNodeID ;
1304             node  [ iRef ] = GMFNode[ aGMFNodeID ];
1305           }
1306         }
1307
1308         switch (token)
1309         {
1310         case GmfEdges:
1311           if (fullyCreatedElement) {
1312             aCreatedElem = theHelper->AddEdge( node[0], node[1], /*id =*/0, /*force3d =*/false );
1313             if (anEdgeGroupByGhs3dId.size() && !anEdgeGroupByGhs3dId[iElem].empty())
1314               addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, anEdgeGroupByGhs3dId[iElem], groupsToRemove);
1315           }
1316           break;
1317         case GmfTriangles:
1318           if (fullyCreatedElement) {
1319             aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], /*id =*/0, /*force3d =*/false );
1320             // for ( int iRef = 0; iRef < nbRef; iRef++ )
1321             //   nodeAssigne[ nodeID[ iRef ]] = 1;
1322             if (aFaceGroupByGhs3dId.size() && !aFaceGroupByGhs3dId[iElem].empty())
1323               addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, aFaceGroupByGhs3dId[iElem], groupsToRemove);
1324           }
1325           break;
1326         case GmfQuadrilaterals:
1327           if (fullyCreatedElement) {
1328             theHelper->AddFace( node[0], node[1], node[2], node[3], /*id =*/0, /*force3d =*/false );
1329             // for ( int iRef = 0; iRef < nbRef; iRef++ )
1330             //   nodeAssigne[ nodeID[ iRef ]] = 1;
1331           }
1332           break;
1333         case GmfTetrahedra:
1334           if ( elemSearcher ) {
1335             // Issue 0020682. Avoid creating nodes and tetras at place where
1336             // volumic elements already exist
1337             if ( !node[1] || !node[0] || !node[2] || !node[3] )
1338               continue;
1339             if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) +
1340                                                     SMESH_TNodeXYZ(node[1]) +
1341                                                     SMESH_TNodeXYZ(node[2]) +
1342                                                     SMESH_TNodeXYZ(node[3]) ) / 4.,
1343                                                     SMDSAbs_Volume, foundVolumes ))
1344             break;
1345           }
1346           theHelper->AddVolume( node[1], node[0], node[2], node[3], /*id =*/0, /*force3d =*/false );
1347 //           theMeshDS->SetMeshElementOnShape( aTet, shapeID );
1348           break;
1349         case GmfHexahedra:
1350           if ( elemSearcher ) {
1351             // Issue 0020682. Avoid creating nodes and tetras at place where
1352             // volumic elements already exist
1353             if ( !node[1] || !node[0] || !node[2] || !node[3] || !node[4] || !node[5] || !node[6] || !node[7])
1354               continue;
1355             if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) +
1356                                                     SMESH_TNodeXYZ(node[1]) +
1357                                                     SMESH_TNodeXYZ(node[2]) +
1358                                                     SMESH_TNodeXYZ(node[3]) +
1359                                                     SMESH_TNodeXYZ(node[4]) +
1360                                                     SMESH_TNodeXYZ(node[5]) +
1361                                                     SMESH_TNodeXYZ(node[6]) +
1362                                                     SMESH_TNodeXYZ(node[7])) / 8.,
1363                                                     SMDSAbs_Volume, foundVolumes ))
1364             break;
1365           }
1366           theHelper->AddVolume( node[0], node[3], node[2], node[1],
1367                                 node[4], node[7], node[6], node[5], /*id =*/0, /*force3d =*/false );
1368 //           theMeshDS->SetMeshElementOnShape( aTet, shapeID );
1369           break;
1370         default: continue;
1371         }
1372       }
1373       break;
1374     }
1375     }
1376   }
1377
1378   // for ( int i = 0; i < nbVertices; ++i ) {
1379   //   if ( !nodeAssigne[ i+1 ])
1380   //     theMeshDS->SetNodeInVolume( GMFNode[ i+1 ], shapeID );
1381   // }
1382
1383   GmfCloseMesh(InpMsh);
1384   delete [] GMFNode;
1385   //delete [] nodeAssigne;
1386 #ifdef _DEBUG_
1387   MESSAGE("Nb subdomains " << subdomainId2tetraId.size());
1388   std::map<int, std::set<int> >::const_iterator subdomainIt = subdomainId2tetraId.begin();
1389   TCollection_AsciiString aSubdomainFileName = theFile;
1390   aSubdomainFileName = aSubdomainFileName + ".subdomain";
1391   ofstream aSubdomainFile  ( aSubdomainFileName.ToCString()  , ios::out);
1392
1393   aSubdomainFile << "Nb subdomains " << subdomainId2tetraId.size() << std::endl;
1394   for(;subdomainIt != subdomainId2tetraId.end() ; ++subdomainIt) {
1395     int subdomainId = subdomainIt->first;
1396     std::set<int> tetraIds = subdomainIt->second;
1397     MESSAGE("Subdomain #"<<subdomainId<<": "<<tetraIds.size()<<" tetrahedrons");
1398     std::set<int>::const_iterator tetraIdsIt = tetraIds.begin();
1399     aSubdomainFile << subdomainId << std::endl;
1400     for(;tetraIdsIt != tetraIds.end() ; ++tetraIdsIt) {
1401       aSubdomainFile << (*tetraIdsIt) << " ";
1402     }
1403     aSubdomainFile << std::endl;
1404   }
1405   aSubdomainFile.close();
1406 #endif
1407
1408   return true;
1409 }
1410
1411 static bool writeGMFFile(const char*                                     theMeshFileName,
1412                          const char*                                     theRequiredFileName,
1413                          const char*                                     theSolFileName,
1414                          const SMESH_ProxyMesh&                          theProxyMesh,
1415                          SMESH_Mesh *                                    theMesh,
1416                          std::vector <const SMDS_MeshNode*> &            theNodeByGhs3dId,
1417                          std::map<const SMDS_MeshNode*,int> &            aNodeToGhs3dIdMap,
1418                          std::vector<std::string> &                      aNodeGroupByGhs3dId,
1419                          std::vector<std::string> &                      anEdgeGroupByGhs3dId,
1420                          std::vector<std::string> &                      aFaceGroupByGhs3dId,
1421                          GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap & theEnforcedNodes,
1422                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges,
1423                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles,
1424                          std::map<std::vector<double>, std::string> &    enfVerticesWithGroup,
1425                          GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices)
1426 {
1427   MESSAGE("writeGMFFile w/o geometry");
1428   std::string tmpStr;
1429   int idx, idxRequired = 0, idxSol = 0;
1430   const int dummyint = 0;
1431   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt;
1432   std::vector<double> enfVertexSizes;
1433   const SMDS_MeshElement* elem;
1434   TIDSortedElemSet anElemSet, theKeptEnforcedEdges, theKeptEnforcedTriangles;
1435   SMDS_ElemIteratorPtr nodeIt;
1436   std::vector <const SMDS_MeshNode*> theEnforcedNodeByGhs3dId;
1437   map<const SMDS_MeshNode*,int> anEnforcedNodeToGhs3dIdMap, anExistingEnforcedNodeToGhs3dIdMap;
1438   std::vector< const SMDS_MeshElement* > foundElems;
1439   map<const SMDS_MeshNode*,TopAbs_State> aNodeToTopAbs_StateMap;
1440   int nbFoundElems;
1441   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap::iterator elemIt;
1442   TIDSortedElemSet::iterator elemSetIt;
1443   bool isOK;
1444   auto_ptr< SMESH_ElementSearcher > pntCls ( SMESH_MeshEditor( theMesh ).GetElementSearcher());
1445
1446   int nbEnforcedVertices = theEnforcedVertices.size();
1447
1448   // count faces
1449   int nbFaces = theProxyMesh.NbFaces();
1450   int nbNodes;
1451
1452   // groups management
1453   int usedEnforcedNodes = 0;
1454   std::string gn = "";
1455
1456   if ( nbFaces == 0 )
1457     return false;
1458
1459   idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1460   if (!idx)
1461     return false;
1462
1463   /* ========================== FACES ========================== */
1464   /* TRIANGLES ========================== */
1465   SMDS_ElemIteratorPtr eIt = theProxyMesh.GetFaces();
1466   while ( eIt->more() )
1467   {
1468     elem = eIt->next();
1469     anElemSet.insert(elem);
1470     nodeIt = elem->nodesIterator();
1471     nbNodes = elem->NbCornerNodes();
1472     while ( nodeIt->more() && nbNodes--)
1473     {
1474       // find GHS3D ID
1475       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1476       int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
1477       aNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1478     }
1479   }
1480
1481   /* EDGES ========================== */
1482
1483   // Iterate over the enforced edges
1484   for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
1485     elem = elemIt->first;
1486     isOK = true;
1487     nodeIt = elem->nodesIterator();
1488     nbNodes = 2;
1489     while ( nodeIt->more() && nbNodes-- ) {
1490       // find GHS3D ID
1491       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1492       // Test if point is inside shape to mesh
1493       gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1494       TopAbs_State result = pntCls->GetPointState( myPoint );
1495       if ( result == TopAbs_OUT ) {
1496         isOK = false;
1497         break;
1498       }
1499       aNodeToTopAbs_StateMap.insert( make_pair( node, result ));
1500     }
1501     if (isOK) {
1502       nodeIt = elem->nodesIterator();
1503       nbNodes = 2;
1504       int newId = -1;
1505       while ( nodeIt->more() && nbNodes-- ) {
1506         // find GHS3D ID
1507         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1508         gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1509         nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1510 #ifdef _DEBUG_
1511         std::cout << "Node at "<<node->X()<<", "<<node->Y()<<", "<<node->Z()<<std::endl;
1512         std::cout << "Nb nodes found : "<<nbFoundElems<<std::endl;
1513 #endif
1514         if (nbFoundElems ==0) {
1515           if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1516             newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
1517             anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1518           }
1519         }
1520         else if (nbFoundElems ==1) {
1521           const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1522           newId = (*aNodeToGhs3dIdMap.find(existingNode)).second;
1523           anExistingEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1524         }
1525         else
1526           isOK = false;
1527 #ifdef _DEBUG_
1528         std::cout << "GHS3D node ID: "<<newId<<std::endl;
1529 #endif
1530       }
1531       if (isOK)
1532         theKeptEnforcedEdges.insert(elem);
1533     }
1534   }
1535
1536   /* ENFORCED TRIANGLES ========================== */
1537
1538   // Iterate over the enforced triangles
1539   for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
1540     elem = elemIt->first;
1541     isOK = true;
1542     nodeIt = elem->nodesIterator();
1543     nbNodes = 3;
1544     while ( nodeIt->more() && nbNodes--) {
1545       // find GHS3D ID
1546       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1547       // Test if point is inside shape to mesh
1548       gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1549       TopAbs_State result = pntCls->GetPointState( myPoint );
1550       if ( result == TopAbs_OUT ) {
1551         isOK = false;
1552         break;
1553       }
1554       aNodeToTopAbs_StateMap.insert( make_pair( node, result ));
1555     }
1556     if (isOK) {
1557       nodeIt = elem->nodesIterator();
1558       nbNodes = 3;
1559       int newId = -1;
1560       while ( nodeIt->more() && nbNodes--) {
1561         // find GHS3D ID
1562         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1563         gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1564         nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1565 #ifdef _DEBUG_
1566         std::cout << "Nb nodes found : "<<nbFoundElems<<std::endl;
1567 #endif
1568         if (nbFoundElems ==0) {
1569           if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1570             newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
1571             anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1572           }
1573         }
1574         else if (nbFoundElems ==1) {
1575           const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1576           newId = (*aNodeToGhs3dIdMap.find(existingNode)).second;
1577           anExistingEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1578         }
1579         else
1580           isOK = false;
1581 #ifdef _DEBUG_
1582         std::cout << "GHS3D node ID: "<<newId<<std::endl;
1583 #endif
1584       }
1585       if (isOK)
1586         theKeptEnforcedTriangles.insert(elem);
1587     }
1588   }
1589
1590   // put nodes to theNodeByGhs3dId vector
1591 #ifdef _DEBUG_
1592   std::cout << "aNodeToGhs3dIdMap.size(): "<<aNodeToGhs3dIdMap.size()<<std::endl;
1593 #endif
1594   theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
1595   map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
1596   for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
1597   {
1598 //     std::cout << "n2id->first: "<<n2id->first<<std::endl;
1599     theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
1600   }
1601
1602   // put nodes to anEnforcedNodeToGhs3dIdMap vector
1603 #ifdef _DEBUG_
1604   std::cout << "anEnforcedNodeToGhs3dIdMap.size(): "<<anEnforcedNodeToGhs3dIdMap.size()<<std::endl;
1605 #endif
1606   theEnforcedNodeByGhs3dId.resize( anEnforcedNodeToGhs3dIdMap.size());
1607   n2id = anEnforcedNodeToGhs3dIdMap.begin();
1608   for ( ; n2id != anEnforcedNodeToGhs3dIdMap.end(); ++ n2id)
1609   {
1610     if (n2id->second > aNodeToGhs3dIdMap.size()) {
1611       theEnforcedNodeByGhs3dId[ n2id->second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // ghs3d ids count from 1
1612     }
1613   }
1614
1615
1616   /* ========================== NODES ========================== */
1617   vector<const SMDS_MeshNode*> theOrderedNodes, theRequiredNodes;
1618   std::set< std::vector<double> > nodesCoords;
1619   vector<const SMDS_MeshNode*>::const_iterator ghs3dNodeIt = theNodeByGhs3dId.begin();
1620   vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
1621
1622   (theNodeByGhs3dId.size() <= 1) ? tmpStr = " node" : " nodes";
1623   std::cout << theNodeByGhs3dId.size() << tmpStr << " from mesh ..." << std::endl;
1624   for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
1625   {
1626     const SMDS_MeshNode* node = *ghs3dNodeIt;
1627     std::vector<double> coords;
1628     coords.push_back(node->X());
1629     coords.push_back(node->Y());
1630     coords.push_back(node->Z());
1631     nodesCoords.insert(coords);
1632     theOrderedNodes.push_back(node);
1633   }
1634
1635   // Iterate over the enforced nodes given by enforced elements
1636   ghs3dNodeIt = theEnforcedNodeByGhs3dId.begin();
1637   after  = theEnforcedNodeByGhs3dId.end();
1638   (theEnforcedNodeByGhs3dId.size() <= 1) ? tmpStr = " node" : " nodes";
1639   std::cout << theEnforcedNodeByGhs3dId.size() << tmpStr << " from enforced elements ..." << std::endl;
1640   for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
1641   {
1642     const SMDS_MeshNode* node = *ghs3dNodeIt;
1643     std::vector<double> coords;
1644     coords.push_back(node->X());
1645     coords.push_back(node->Y());
1646     coords.push_back(node->Z());
1647 #ifdef _DEBUG_
1648     std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z();
1649 #endif
1650
1651     if (nodesCoords.find(coords) != nodesCoords.end()) {
1652       // node already exists in original mesh
1653 #ifdef _DEBUG_
1654       std::cout << " found" << std::endl;
1655 #endif
1656       continue;
1657     }
1658
1659     if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) {
1660       // node already exists in enforced vertices
1661 #ifdef _DEBUG_
1662       std::cout << " found" << std::endl;
1663 #endif
1664       continue;
1665     }
1666
1667 //     gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1668 //     nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1669 //     if (nbFoundElems ==0) {
1670 //       std::cout << " not found" << std::endl;
1671 //       if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1672 //         nodesCoords.insert(coords);
1673 //         theOrderedNodes.push_back(node);
1674 //       }
1675 //     }
1676 //     else {
1677 //       std::cout << " found in initial mesh" << std::endl;
1678 //       const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1679 //       nodesCoords.insert(coords);
1680 //       theOrderedNodes.push_back(existingNode);
1681 //     }
1682
1683 #ifdef _DEBUG_
1684     std::cout << " not found" << std::endl;
1685 #endif
1686
1687     nodesCoords.insert(coords);
1688     theOrderedNodes.push_back(node);
1689 //     theRequiredNodes.push_back(node);
1690   }
1691
1692
1693   // Iterate over the enforced nodes
1694   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator enfNodeIt;
1695   (theEnforcedNodes.size() <= 1) ? tmpStr = " node" : " nodes";
1696   std::cout << theEnforcedNodes.size() << tmpStr << " from enforced nodes ..." << std::endl;
1697   for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt)
1698   {
1699     const SMDS_MeshNode* node = enfNodeIt->first;
1700     std::vector<double> coords;
1701     coords.push_back(node->X());
1702     coords.push_back(node->Y());
1703     coords.push_back(node->Z());
1704 #ifdef _DEBUG_
1705     std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z();
1706 #endif
1707
1708     // Test if point is inside shape to mesh
1709     gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1710     TopAbs_State result = pntCls->GetPointState( myPoint );
1711     if ( result == TopAbs_OUT ) {
1712 #ifdef _DEBUG_
1713       std::cout << " out of volume" << std::endl;
1714 #endif
1715       continue;
1716     }
1717
1718     if (nodesCoords.find(coords) != nodesCoords.end()) {
1719 #ifdef _DEBUG_
1720       std::cout << " found in nodesCoords" << std::endl;
1721 #endif
1722 //       theRequiredNodes.push_back(node);
1723       continue;
1724     }
1725
1726     if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) {
1727 #ifdef _DEBUG_
1728       std::cout << " found in theEnforcedVertices" << std::endl;
1729 #endif
1730       continue;
1731     }
1732
1733 //     nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1734 //     if (nbFoundElems ==0) {
1735 //       std::cout << " not found" << std::endl;
1736 //       if (result == TopAbs_IN) {
1737 //         nodesCoords.insert(coords);
1738 //         theRequiredNodes.push_back(node);
1739 //       }
1740 //     }
1741 //     else {
1742 //       std::cout << " found in initial mesh" << std::endl;
1743 //       const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1744 // //       nodesCoords.insert(coords);
1745 //       theRequiredNodes.push_back(existingNode);
1746 //     }
1747 //
1748 //
1749 //
1750 //     if (pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems) == 0)
1751 //       continue;
1752
1753 //     if ( result != TopAbs_IN )
1754 //       continue;
1755
1756 #ifdef _DEBUG_
1757     std::cout << " not found" << std::endl;
1758 #endif
1759     nodesCoords.insert(coords);
1760 //     theOrderedNodes.push_back(node);
1761     theRequiredNodes.push_back(node);
1762   }
1763   int requiredNodes = theRequiredNodes.size();
1764
1765   int solSize = 0;
1766   std::vector<std::vector<double> > ReqVerTab;
1767   if (nbEnforcedVertices) {
1768 //    ReqVerTab.clear();
1769     (nbEnforcedVertices <= 1) ? tmpStr = " node" : " nodes";
1770     std::cout << nbEnforcedVertices << tmpStr << " from enforced vertices ..." << std::endl;
1771     // Iterate over the enforced vertices
1772     for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
1773       double x = vertexIt->first[0];
1774       double y = vertexIt->first[1];
1775       double z = vertexIt->first[2];
1776       // Test if point is inside shape to mesh
1777       gp_Pnt myPoint(x,y,z);
1778       TopAbs_State result = pntCls->GetPointState( myPoint );
1779       if ( result == TopAbs_OUT )
1780         continue;
1781       //if (pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems) == 0)
1782       //continue;
1783
1784 //       if ( result != TopAbs_IN )
1785 //         continue;
1786       std::vector<double> coords;
1787       coords.push_back(x);
1788       coords.push_back(y);
1789       coords.push_back(z);
1790       ReqVerTab.push_back(coords);
1791       enfVertexSizes.push_back(vertexIt->second);
1792       solSize++;
1793     }
1794   }
1795
1796
1797   // GmfVertices
1798   std::cout << "Begin writting required nodes in GmfVertices" << std::endl;
1799   std::cout << "Nb vertices: " << theOrderedNodes.size() << std::endl;
1800   GmfSetKwd(idx, GmfVertices, theOrderedNodes.size()/*+solSize*/);
1801   for (ghs3dNodeIt = theOrderedNodes.begin();ghs3dNodeIt != theOrderedNodes.end();++ghs3dNodeIt) {
1802     GmfSetLin(idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint);
1803   }
1804
1805   std::cout << "End writting required nodes in GmfVertices" << std::endl;
1806
1807   if (requiredNodes + solSize) {
1808     std::cout << "Begin writting in req and sol file" << std::endl;
1809     aNodeGroupByGhs3dId.resize( requiredNodes + solSize );
1810     idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1811     if (!idxRequired) {
1812       GmfCloseMesh(idx);
1813       return false;
1814     }
1815     idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1816     if (!idxSol) {
1817       GmfCloseMesh(idx);
1818       if (idxRequired)
1819         GmfCloseMesh(idxRequired);
1820       return false;
1821     }
1822     int TypTab[] = {GmfSca};
1823     double ValTab[] = {0.0};
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, ValTab);
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       << -1 << 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 }
4315
4316 namespace
4317 {
4318   //================================================================================
4319   /*!
4320    * \brief Sub-mesh event listener setting enforced elements as soon as an enforced
4321    *        mesh is loaded
4322    */
4323   struct _EnforcedMeshRestorer : public SMESH_subMeshEventListener
4324   {
4325     _EnforcedMeshRestorer():
4326       SMESH_subMeshEventListener( /*isDeletable = */true, Name() )
4327     {}
4328
4329     //================================================================================
4330     /*!
4331      * \brief Returns an ID of listener
4332      */
4333     static const char* Name() { return "GHS3DPlugin_GHS3D::_EnforcedMeshRestorer"; }
4334
4335     //================================================================================
4336     /*!
4337      * \brief Treat events of the subMesh
4338      */
4339     void ProcessEvent(const int                       event,
4340                       const int                       eventType,
4341                       SMESH_subMesh*                  subMesh,
4342                       SMESH_subMeshEventListenerData* data,
4343                       const SMESH_Hypothesis*         hyp)
4344     {
4345       if ( SMESH_subMesh::SUBMESH_LOADED == event &&
4346            SMESH_subMesh::COMPUTE_EVENT  == eventType &&
4347            data &&
4348            !data->mySubMeshes.empty() )
4349       {
4350         // An enforced mesh (subMesh->_father) has been loaded from hdf file
4351         if ( GHS3DPlugin_Hypothesis* hyp = GetGHSHypothesis( data->mySubMeshes.front() ))
4352           hyp->RestoreEnfElemsByMeshes();
4353       }
4354     }
4355     //================================================================================
4356     /*!
4357      * \brief Returns GHS3DPlugin_Hypothesis used to compute a subMesh
4358      */
4359     static GHS3DPlugin_Hypothesis* GetGHSHypothesis( SMESH_subMesh* subMesh )
4360     {
4361       SMESH_HypoFilter ghsHypFilter( SMESH_HypoFilter::HasName( "GHS3D_Parameters" ));
4362       return (GHS3DPlugin_Hypothesis* )
4363         subMesh->GetFather()->GetHypothesis( subMesh->GetSubShape(),
4364                                              ghsHypFilter,
4365                                              /*visitAncestors=*/true);
4366     }
4367   };
4368 }
4369
4370 //================================================================================
4371 /*!
4372  * \brief Set an event listener to set enforced elements as soon as an enforced
4373  *        mesh is loaded
4374  */
4375 //================================================================================
4376
4377 void GHS3DPlugin_GHS3D::SubmeshRestored(SMESH_subMesh* subMesh)
4378 {
4379   if ( GHS3DPlugin_Hypothesis* hyp = _EnforcedMeshRestorer::GetGHSHypothesis( subMesh ))
4380   {
4381     GHS3DPlugin_Hypothesis::TGHS3DEnforcedMeshList enfMeshes = hyp->_GetEnforcedMeshes();
4382     GHS3DPlugin_Hypothesis::TGHS3DEnforcedMeshList::iterator it = enfMeshes.begin();
4383     for(;it != enfMeshes.end();++it) {
4384       GHS3DPlugin_Hypothesis::TGHS3DEnforcedMesh* enfMesh = *it;
4385       if ( SMESH_Mesh* mesh = GetMeshByPersistentID( enfMesh->persistID ))
4386       {
4387         SMESH_subMesh* smToListen = mesh->GetSubMesh( mesh->GetShapeToMesh() );
4388         // a listener set to smToListen will care of hypothesis stored in SMESH_EventListenerData
4389         subMesh->SetEventListener( new _EnforcedMeshRestorer(),
4390                                    SMESH_subMeshEventListenerData::MakeData( subMesh ),
4391                                    smToListen);
4392       }
4393     }
4394   }
4395 }