src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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