src/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
  31 #include <Basics_Utils.hxx>
  32
  33 #include "SMESH_Gen.hxx"
  34 #include "SMESH_Mesh.hxx"
  35 #include "SMESH_Comment.hxx"
  36 #include "SMESH_MesherHelper.hxx"
  37 #include "SMESH_MeshEditor.hxx"
  38
  39 #include "SMDS_MeshElement.hxx"
  40 #include "SMDS_MeshNode.hxx"
  41 #include "SMDS_FaceOfNodes.hxx"
  42 #include "SMDS_VolumeOfNodes.hxx"
  43
  44 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  45 #include <StdMeshers_ViscousLayers.hxx>
  46
  47 #include <BRepAdaptor_Surface.hxx>
  48 #include <BRepBndLib.hxx>
  49 #include <BRepClass3d_SolidClassifier.hxx>
  50 #include <BRepTools.hxx>
  51 #include <BRep_Tool.hxx>
  52 #include <Bnd_Box.hxx>
  53 #include <GeomAPI_ProjectPointOnSurf.hxx>
  54 #include <OSD_File.hxx>
  55 #include <Precision.hxx>
  56 #include <Quantity_Parameter.hxx>
  57 #include <Standard_ProgramError.hxx>
  58 #include <Standard_ErrorHandler.hxx>
  59 #include <Standard_Failure.hxx>
  60 #include <TopExp.hxx>
  61 #include <TopExp_Explorer.hxx>
  62 #include <TopTools_IndexedMapOfShape.hxx>
  63 #include <TopTools_ListIteratorOfListOfShape.hxx>
  64 #include <TopoDS.hxx>
  65 //#include <BRepClass_FaceClassifier.hxx>
  66 #include <TopTools_MapOfShape.hxx>
  67 #include <BRepGProp.hxx>
  68 #include <GProp_GProps.hxx>
  69
  70 #include "utilities.h"
  71
  72 #ifdef WIN32
  73 #include <io.h>
  74 #else
  75 #include <sys/sysinfo.h>
  76 #endif
  77
  78 using namespace std;
  79
  80 //#include <Standard_Stream.hxx>
  81
  82
  83 #define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
  84
  85 #ifdef _DEBUG_
  86 #define DUMP(txt) \
  87 //  std::cout << txt
  88 #else
  89 #define DUMP(txt)
  90 #endif
  91
  92 extern "C"
  93 {
  94 #ifndef WNT
  95 #include <unistd.h>
  96 #include <sys/mman.h>
  97 #endif
  98 #include <sys/stat.h>
  99 #include <fcntl.h>
 100 }
 101
 102 #define HOLE_ID -1
 103
 104 typedef const list<const SMDS_MeshFace*> TTriaList;
 105
 106 static void removeFile( const TCollection_AsciiString& fileName )
 107 {
 108   try {
 109     OSD_File( fileName ).Remove();
 110   }
 111   catch ( Standard_ProgramError ) {
 112     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
 113   }
 114 }
 115
 116 //=============================================================================
 117 /*!
 118  *
 119  */
 120 //=============================================================================
 121
 122 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
 123   : SMESH_3D_Algo(hypId, studyId, gen)
 124 {
 125   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
 126   _name = "GHS3D_3D";
 127   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
 128   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
 129   _iShape=0;
 130   _nbShape=0;
 131   _compatibleHypothesis.push_back("GHS3D_Parameters");
 132   _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() );
 133   _requireShape = false; // can work without shape
 134 }
 135
 136 //=============================================================================
 137 /*!
 138  *
 139  */
 140 //=============================================================================
 141
 142 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
 143 {
 144   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
 145 }
 146
 147 //=============================================================================
 148 /*!
 149  *
 150  */
 151 //=============================================================================
 152
 153 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&         aMesh,
 154                                           const TopoDS_Shape& aShape,
 155                                           Hypothesis_Status&  aStatus )
 156 {
 157   aStatus = SMESH_Hypothesis::HYP_OK;
 158
 159   _hyp = 0;
 160   _viscousLayersHyp = 0;
 161   _keepFiles = false;
 162
 163   const list <const SMESHDS_Hypothesis * >& hyps =
 164     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 165   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 166   for ( ; h != hyps.end(); ++h )
 167   {
 168     if ( !_hyp )
 169       _hyp = dynamic_cast< const GHS3DPlugin_Hypothesis*> ( *h );
 170     if ( !_viscousLayersHyp )
 171       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 172   }
 173   if ( _hyp )
 174     _keepFiles = _hyp->GetKeepFiles();
 175
 176   return true;
 177 }
 178
 179 //=======================================================================
 180 //function : findShape
 181 //purpose  :
 182 //=======================================================================
 183
 184 static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[],
 185                               TopoDS_Shape        aShape,
 186                               const TopoDS_Shape  shape[],
 187                               double**            box,
 188                               const int           nShape,
 189                               TopAbs_State *      state = 0)
 190 {
 191   gp_XYZ aPnt(0,0,0);
 192   int j, iShape, nbNode = 4;
 193
 194   for ( j=0; j<nbNode; j++ ) {
 195     gp_XYZ p ( aNode[j]->X(), aNode[j]->Y(), aNode[j]->Z() );
 196     if ( aNode[j]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) {
 197       aPnt = p;
 198       break;
 199     }
 200     aPnt += p / nbNode;
 201   }
 202
 203   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 204   if (state) *state = SC.State();
 205   if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) {
 206     for (iShape = 0; iShape < nShape; iShape++) {
 207       aShape = shape[iShape];
 208       if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() ||
 209               aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() ||
 210               aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) {
 211         BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 212         if (state) *state = SC.State();
 213         if (SC.State() == TopAbs_IN)
 214           break;
 215       }
 216     }
 217   }
 218   return aShape;
 219 }
 220
 221 //=======================================================================
 222 //function : readMapIntLine
 223 //purpose  :
 224 //=======================================================================
 225
 226 static char* readMapIntLine(char* ptr, int tab[]) {
 227   long int intVal;
 228   std::cout << std::endl;
 229
 230   for ( int i=0; i<17; i++ ) {
 231     intVal = strtol(ptr, &ptr, 10);
 232     if ( i < 3 )
 233       tab[i] = intVal;
 234   }
 235   return ptr;
 236 }
 237
 238 //=======================================================================
 239 //function : countShape
 240 //purpose  :
 241 //=======================================================================
 242
 243 // template < class Mesh, class Shape >
 244 // static int countShape( Mesh* mesh, Shape shape ) {
 245 //   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
 246 //   int nbShape = 0;
 247 //   for ( ; expShape.More(); expShape.Next() ) {
 248 //       nbShape++;
 249 //   }
 250 //   return nbShape;
 251 // }
 252
 253 //=======================================================================
 254 //function : writeFaces
 255 //purpose  :
 256 //=======================================================================
 257
 258 static bool writeFaces (ofstream &             theFile,
 259                         const SMESH_ProxyMesh& theMesh,
 260                         const TopoDS_Shape&    theShape,
 261                         const map <int,int> &  theSmdsToGhs3dIdMap)
 262 {
 263   // record structure:
 264   //
 265   // NB_ELEMS DUMMY_INT
 266   // Loop from 1 to NB_ELEMS
 267   // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
 268
 269   TopoDS_Shape aShape;
 270   const SMESHDS_SubMesh* theSubMesh;
 271   const SMDS_MeshElement* aFace;
 272   const char* space    = "  ";
 273   const int   dummyint = 0;
 274   map<int,int>::const_iterator itOnMap;
 275   SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
 276   int nbNodes, aSmdsID;
 277
 278   // count triangles bound to geometry
 279   int nbTriangles = 0;
 280
 281   TopTools_IndexedMapOfShape facesMap;
 282   TopExp::MapShapes( theShape, TopAbs_FACE, facesMap );
 283
 284   for ( int i = 1; i <= facesMap.Extent(); ++i )
 285     if (( theSubMesh  = theMesh.GetSubMesh( facesMap(i))))
 286       nbTriangles += theSubMesh->NbElements();
 287
 288   std::cout << "    " << facesMap.Extent() << " shapes of 2D dimension" << std::endl;
 289   std::cout << std::endl;
 290
 291   theFile << space << nbTriangles << space << dummyint << std::endl;
 292
 293   for ( int i = 1; i <= facesMap.Extent(); i++ )
 294   {
 295     aShape = facesMap(i);
 296     theSubMesh = theMesh.GetSubMesh(aShape);
 297     if ( !theSubMesh ) continue;
 298     itOnSubMesh = theSubMesh->GetElements();
 299     while ( itOnSubMesh->more() )
 300     {
 301       aFace = itOnSubMesh->next();
 302       nbNodes = aFace->NbNodes();
 303
 304       theFile << space << nbNodes;
 305
 306       itOnSubFace = aFace->nodesIterator();
 307       while ( itOnSubFace->more() ) {
 308         // find GHS3D ID
 309         aSmdsID = itOnSubFace->next()->GetID();
 310         itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
 311         // if ( itOnMap == theSmdsToGhs3dIdMap.end() ) {
 312         //   cout << "not found node: " << aSmdsID << endl;
 313         //   return false;
 314         // }
 315         ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
 316
 317         theFile << space << (*itOnMap).second;
 318       }
 319
 320       // (NB_NODES + 1) times: DUMMY_INT
 321       for ( int j=0; j<=nbNodes; j++)
 322         theFile << space << dummyint;
 323
 324       theFile << std::endl;
 325     }
 326   }
 327
 328   return true;
 329 }
 330
 331 //=======================================================================
 332 //function : writeFaces
 333 //purpose  : Write Faces in case if generate 3D mesh w/o geometry
 334 //=======================================================================
 335
 336 static bool writeFaces (ofstream &                      theFile,
 337                         const SMESH_ProxyMesh&          theMesh,
 338                         vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
 339 {
 340   // record structure:
 341   //
 342   // NB_ELEMS DUMMY_INT
 343   // Loop from 1 to NB_ELEMS
 344   //   NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
 345
 346   int nbNodes, nbTriangles = 0;
 347   map< const SMDS_MeshNode*,int >::iterator it;
 348   SMDS_ElemIteratorPtr nodeIt;
 349   const SMDS_MeshElement* elem;
 350
 351   const char* space    = "  ";
 352   const int   dummyint = 0;
 353
 354   // count faces
 355   nbTriangles = theMesh.NbFaces();
 356
 357   if ( nbTriangles == 0 )
 358     return false;
 359
 360   std::cout << "The initial 2D mesh contains " << nbTriangles << " faces and ";
 361
 362   // NB_ELEMS DUMMY_INT
 363   theFile << space << nbTriangles << space << dummyint << std::endl;
 364
 365
 366   map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
 367
 368   // Loop from 1 to NB_ELEMS
 369
 370   SMDS_ElemIteratorPtr eIt = theMesh.GetFaces();
 371   while ( eIt->more() )
 372   {
 373     elem = eIt->next();
 374     // NB_NODES PER FACE
 375     nbNodes = elem->NbNodes();
 376     theFile << space << nbNodes;
 377
 378     // NODE_NB_1 NODE_NB_2 ...
 379     nodeIt = elem->nodesIterator();
 380     while ( nodeIt->more() )
 381     {
 382       // find GHS3D ID
 383       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
 384       int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
 385       it = aNodeToGhs3dIdMap.insert( make_pair( node, newId )).first;
 386       theFile << space << it->second;
 387     }
 388
 389     // (NB_NODES + 1) times: DUMMY_INT
 390     for ( int i=0; i<=nbNodes; i++)
 391       theFile << space << dummyint;
 392     theFile << std::endl;
 393   }
 394
 395   // put nodes to theNodeByGhs3dId vector
 396   theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
 397   map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
 398   for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
 399   {
 400     theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
 401   }
 402
 403   return true;
 404 }
 405
 406 //=======================================================================
 407 //function : writePoints
 408 //purpose  :
 409 //=======================================================================
 410
 411 static bool writePoints (ofstream &                       theFile,
 412                          SMESH_MesherHelper&              theHelper,
 413                          map <int,int> &                  theSmdsToGhs3dIdMap,
 414                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap,
 415                          map<vector<double>,double> &     theEnforcedVertices)
 416 {
 417   // record structure:
 418   //
 419   // NB_NODES
 420   // Loop from 1 to NB_NODES
 421   //   X Y Z DUMMY_INT
 422
 423   SMESHDS_Mesh * theMesh = theHelper.GetMeshDS();
 424   int nbNodes = theMesh->NbNodes();
 425   if ( nbNodes == 0 )
 426     return false;
 427   int nbEnforcedVertices = theEnforcedVertices.size();
 428
 429   // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
 430   // The problem is in nodes on degenerated edges, we need to skip them
 431   if ( theHelper.HasDegeneratedEdges() )
 432   {
 433     // here we decrease total nb of nodes by nb of nodes on degenerated edges
 434     set<int> checkedSM;
 435     for (TopExp_Explorer e(theMesh->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
 436     {
 437       SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
 438       if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() )) {
 439         if ( sm->GetSubMeshDS() )
 440           nbNodes -= sm->GetSubMeshDS()->NbNodes();
 441       }
 442     }
 443   }
 444   const char* space    = "  ";
 445   const int   dummyint = 0;
 446
 447   int aGhs3dID = 1;
 448   SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
 449   const SMDS_MeshNode* node;
 450
 451   // NB_NODES
 452   std::cout << std::endl;
 453   std::cout << "The initial 2D mesh contains :" << std::endl;
 454   std::cout << "    " << nbNodes << " nodes" << std::endl;
 455   if (nbEnforcedVertices > 0) {
 456     std::cout << "    " << nbEnforcedVertices << " enforced vertices" << std::endl;
 457   }
 458   std::cout << std::endl;
 459   std::cout << "Start writing in 'points' file ..." << std::endl;
 460   theFile << space << nbNodes << std::endl;
 461
 462   // Loop from 1 to NB_NODES
 463
 464   while ( it->more() )
 465   {
 466     node = it->next();
 467     if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 468          theHelper.IsDegenShape( node->getshapeId() )) // Issue 020674
 469       continue;
 470
 471     theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
 472     theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
 473     aGhs3dID++;
 474
 475     // X Y Z DUMMY_INT
 476     theFile
 477     << space << node->X()
 478     << space << node->Y()
 479     << space << node->Z()
 480     << space << dummyint;
 481
 482     theFile << std::endl;
 483
 484   }
 485
 486   // Iterate over the enforced vertices
 487   GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
 488   const TopoDS_Shape shapeToMesh = theMesh->ShapeToMesh();
 489   for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
 490     double x = vertexIt->first[0];
 491     double y = vertexIt->first[1];
 492     double z = vertexIt->first[2];
 493     // Test if point is inside shape to mesh
 494     gp_Pnt myPoint(x,y,z);
 495     BRepClass3d_SolidClassifier scl(shapeToMesh);
 496     scl.Perform(myPoint, 1e-7);
 497     TopAbs_State result = scl.State();
 498     if ( result == TopAbs_IN ) {
 499         MESSAGE("Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second);
 500         // X Y Z PHY_SIZE DUMMY_INT
 501         theFile
 502         << space << x
 503         << space << y
 504         << space << z
 505         << space << vertexIt->second
 506         << space << dummyint;
 507
 508         theFile << std::endl;
 509     }
 510     else
 511         MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added ");
 512   }
 513
 514
 515   std::cout << std::endl;
 516   std::cout << "End writing in 'points' file." << std::endl;
 517
 518   return true;
 519 }
 520
 521 //=======================================================================
 522 //function : writePoints
 523 //purpose  :
 524 //=======================================================================
 525
 526 static bool writePoints (ofstream &                            theFile,
 527                          SMESH_Mesh *                          theMesh,
 528                          const vector <const SMDS_MeshNode*> & theNodeByGhs3dId,
 529                          map<vector<double>,double> &          theEnforcedVertices)
 530 {
 531   // record structure:
 532   //
 533   // NB_NODES
 534   // Loop from 1 to NB_NODES
 535   //   X Y Z DUMMY_INT
 536
 537   int nbNodes = theNodeByGhs3dId.size();
 538   if ( nbNodes == 0 )
 539     return false;
 540
 541   int nbEnforcedVertices = theEnforcedVertices.size();
 542
 543   const char* space    = "  ";
 544   const int   dummyint = 0;
 545
 546   const SMDS_MeshNode* node;
 547
 548   // NB_NODES
 549   std::cout << std::endl;
 550   std::cout << "The initial 2D mesh contains :" << std::endl;
 551   std::cout << "    " << nbNodes << " nodes" << std::endl;
 552   std::cout << "    " << nbEnforcedVertices << " enforced vertices" << std::endl;
 553   std::cout << std::endl;
 554   std::cout << "Start writing in 'points' file ..." << std::endl;
 555   theFile << space << nbNodes << std::endl;
 556
 557   // Loop from 1 to NB_NODES
 558
 559   vector<const SMDS_MeshNode*>::const_iterator nodeIt = theNodeByGhs3dId.begin();
 560   vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
 561   for ( ; nodeIt != after; ++nodeIt )
 562   {
 563     node = *nodeIt;
 564
 565     // X Y Z DUMMY_INT
 566     theFile
 567     << space << node->X()
 568     << space << node->Y()
 569     << space << node->Z()
 570     << space << dummyint;
 571
 572     theFile << std::endl;
 573
 574   }
 575
 576   // Iterate over the enforced vertices
 577   GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
 578   auto_ptr< SMESH_ElementSearcher > pntCls ( SMESH_MeshEditor( theMesh ).GetElementSearcher());
 579   for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
 580     double x = vertexIt->first[0];
 581     double y = vertexIt->first[1];
 582     double z = vertexIt->first[2];
 583     // Test if point is inside shape to mesh
 584     gp_Pnt myPoint(x,y,z);
 585     TopAbs_State result = pntCls->GetPointState( myPoint );
 586     if ( result == TopAbs_IN ) {
 587         std::cout << "Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second << std::endl;
 588
 589         // X Y Z PHY_SIZE DUMMY_INT
 590         theFile
 591         << space << x
 592         << space << y
 593         << space << z
 594         << space << vertexIt->second
 595         << space << dummyint;
 596
 597         theFile << std::endl;
 598     }
 599   }
 600   std::cout << std::endl;
 601   std::cout << "End writing in 'points' file." << std::endl;
 602
 603   return true;
 604 }
 605
 606 //================================================================================
 607 /*!
 608  * \brief returns true if a triangle defined by the nodes is a temporary face on a
 609  * side facet of pyramid and defines sub-domian inside the pyramid
 610  */
 611 //================================================================================
 612
 613 static bool isTmpFace(const SMDS_MeshNode* node1,
 614                       const SMDS_MeshNode* node2,
 615                       const SMDS_MeshNode* node3)
 616 {
 617   // find a pyramid sharing the 3 nodes
 618   //const SMDS_MeshElement* pyram = 0;
 619   SMDS_ElemIteratorPtr vIt1 = node1->GetInverseElementIterator(SMDSAbs_Volume);
 620   while ( vIt1->more() )
 621   {
 622     const SMDS_MeshElement* pyram = vIt1->next();
 623     if ( pyram->NbCornerNodes() != 5 ) continue;
 624     int i2, i3;
 625     if ( (i2 = pyram->GetNodeIndex( node2 )) >= 0 &&
 626          (i3 = pyram->GetNodeIndex( node3 )) >= 0 )
 627     {
 628       // Triangle defines sub-domian inside the pyramid if it's
 629       // normal points out of the pyram
 630
 631       // make i2 and i3 hold indices of base nodes of the pyram while
 632       // keeping the nodes order in the triangle
 633       const int iApex = 4;
 634       if ( i2 == iApex )
 635         i2 = i3, i3 = pyram->GetNodeIndex( node1 );
 636       else if ( i3 == iApex )
 637         i3 = i2, i2 = pyram->GetNodeIndex( node1 );
 638
 639       int i3base = (i2+1) % 4; // next index after i2 within the pyramid base
 640       return ( i3base != i3 );
 641     }
 642   }
 643   return false;
 644 }
 645
 646 //=======================================================================
 647 //function : findShapeID
 648 //purpose  : find the solid corresponding to GHS3D sub-domain following
 649 //           the technique proposed in GHS3D manual (available within
 650 //           ghs3d installation) in chapter "B.4 Subdomain (sub-region) assignment".
 651 //           In brief: normal of the triangle defined by the given nodes
 652 //           points out of the domain it is associated to
 653 //=======================================================================
 654
 655 static int findShapeID(SMESH_Mesh&          mesh,
 656                        const SMDS_MeshNode* node1,
 657                        const SMDS_MeshNode* node2,
 658                        const SMDS_MeshNode* node3,
 659                        const bool           toMeshHoles)
 660 {
 661   const int invalidID = 0;
 662   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 663
 664   // face the nodes belong to
 665   const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
 666   if ( !face )
 667     return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
 668 #ifdef _DEBUG_
 669   std::cout << "bnd face " << face->GetID() << " - ";
 670 #endif
 671   // geom face the face assigned to
 672   SMESH_MeshEditor editor(&mesh);
 673   int geomFaceID = editor.FindShape( face );
 674   if ( !geomFaceID )
 675     return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
 676   TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
 677   if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
 678     return invalidID;
 679   TopoDS_Face geomFace = TopoDS::Face( shape );
 680
 681   // solids bounded by geom face
 682   TopTools_IndexedMapOfShape solids, shells;
 683   TopTools_ListIteratorOfListOfShape ansIt = mesh.GetAncestors(geomFace);
 684   for ( ; ansIt.More(); ansIt.Next() ) {
 685     switch ( ansIt.Value().ShapeType() ) {
 686     case TopAbs_SOLID:
 687       solids.Add( ansIt.Value() ); break;
 688     case TopAbs_SHELL:
 689       shells.Add( ansIt.Value() ); break;
 690     default:;
 691     }
 692   }
 693   // analyse found solids
 694   if ( solids.Extent() == 0 || shells.Extent() == 0)
 695     return invalidID;
 696
 697   const TopoDS_Solid& solid1 = TopoDS::Solid( solids(1) );
 698   if ( solids.Extent() == 1 )
 699   {
 700     if ( toMeshHoles )
 701       return meshDS->ShapeToIndex( solid1 );
 702
 703     // - Are we at a hole boundary face?
 704     if ( shells(1).IsSame( BRepTools::OuterShell( solid1 )) )
 705     { // - No, but maybe a hole is bound by two shapes? Does shells(1) touches another shell?
 706       bool touch = false;
 707       TopExp_Explorer eExp( shells(1), TopAbs_EDGE );
 708       // check if any edge of shells(1) belongs to another shell
 709       for ( ; eExp.More() && !touch; eExp.Next() ) {
 710         ansIt = mesh.GetAncestors( eExp.Current() );
 711         for ( ; ansIt.More() && !touch; ansIt.Next() ) {
 712           if ( ansIt.Value().ShapeType() == TopAbs_SHELL )
 713             touch = ( !ansIt.Value().IsSame( shells(1) ));
 714         }
 715       }
 716       if (!touch)
 717         return meshDS->ShapeToIndex( solid1 );
 718     }
 719   }
 720   // find orientation of geom face within the first solid
 721   TopExp_Explorer fExp( solid1, TopAbs_FACE );
 722   for ( ; fExp.More(); fExp.Next() )
 723     if ( geomFace.IsSame( fExp.Current() )) {
 724       geomFace = TopoDS::Face( fExp.Current() );
 725       break;
 726     }
 727   if ( !fExp.More() )
 728     return invalidID; // face not found
 729
 730   // normale to triangle
 731   gp_Pnt node1Pnt ( node1->X(), node1->Y(), node1->Z() );
 732   gp_Pnt node2Pnt ( node2->X(), node2->Y(), node2->Z() );
 733   gp_Pnt node3Pnt ( node3->X(), node3->Y(), node3->Z() );
 734   gp_Vec vec12( node1Pnt, node2Pnt );
 735   gp_Vec vec13( node1Pnt, node3Pnt );
 736   gp_Vec meshNormal = vec12 ^ vec13;
 737   if ( meshNormal.SquareMagnitude() < DBL_MIN )
 738     return invalidID;
 739
 740   // get normale to geomFace at any node
 741   bool geomNormalOK = false;
 742   gp_Vec geomNormal;
 743   const SMDS_MeshNode* nodes[3] = { node1, node2, node3 };
 744   SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace );
 745   for ( int i = 0; !geomNormalOK && i < 3; ++i )
 746   {
 747     // find UV of i-th node on geomFace
 748     const SMDS_MeshNode* nNotOnSeamEdge = 0;
 749     if ( helper.IsSeamShape( nodes[i]->getshapeId() ))
 750       if ( helper.IsSeamShape( nodes[(i+1)%3]->getshapeId() ))
 751         nNotOnSeamEdge = nodes[(i+2)%3];
 752       else
 753         nNotOnSeamEdge = nodes[(i+1)%3];
 754     bool uvOK;
 755     gp_XY uv = helper.GetNodeUV( geomFace, nodes[i], nNotOnSeamEdge, &uvOK );
 756     // check that uv is correct
 757     if (uvOK) {
 758       double tol = 1e-6;
 759       TopoDS_Shape nodeShape = helper.GetSubShapeByNode( nodes[i], meshDS );
 760       if ( !nodeShape.IsNull() )
 761         switch ( nodeShape.ShapeType() )
 762         {
 763         case TopAbs_FACE:   tol = BRep_Tool::Tolerance( TopoDS::Face( nodeShape )); break;
 764         case TopAbs_EDGE:   tol = BRep_Tool::Tolerance( TopoDS::Edge( nodeShape )); break;
 765         case TopAbs_VERTEX: tol = BRep_Tool::Tolerance( TopoDS::Vertex( nodeShape )); break;
 766         default:;
 767         }
 768       gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
 769       BRepAdaptor_Surface surface( geomFace );
 770       uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
 771       if ( uvOK ) {
 772         // normale to geomFace at UV
 773         gp_Vec du, dv;
 774         surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
 775         geomNormal = du ^ dv;
 776         if ( geomFace.Orientation() == TopAbs_REVERSED )
 777           geomNormal.Reverse();
 778         geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
 779       }
 780     }
 781   }
 782   if ( !geomNormalOK)
 783     return invalidID;
 784
 785   // compare normals
 786   bool isReverse = ( meshNormal * geomNormal ) < 0;
 787   if ( !isReverse )
 788     return meshDS->ShapeToIndex( solid1 );
 789
 790   if ( solids.Extent() == 1 )
 791     return HOLE_ID; // we are inside a hole
 792   else
 793     return meshDS->ShapeToIndex( solids(2) );
 794 }
 795
 796 //=======================================================================
 797 //function : readResultFile
 798 //purpose  :
 799 //=======================================================================
 800
 801 static bool readResultFile(const int                       fileOpen,
 802 #ifdef WNT
 803                            const char*                     fileName,
 804 #endif
 805                            SMESH_Mesh&                     theMesh,
 806                            TopoDS_Shape                    tabShape[],
 807                            double**                        tabBox,
 808                            const int                       nbShape,
 809                            map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
 810                            bool                            toMeshHoles,
 811                            int                             nbEnforcedVertices)
 812 {
 813   MESSAGE("GHS3DPlugin_GHS3D::readResultFile()");
 814   Kernel_Utils::Localizer loc;
 815   struct stat status;
 816   size_t      length;
 817
 818   char *ptr, *mapPtr;
 819   char *tetraPtr;
 820   char *shapePtr;
 821
 822   SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
 823
 824   int fileStat;
 825   int nbElems, nbNodes, nbInputNodes;
 826   int nodeId/*, triangleId*/;
 827   int nbTriangle;
 828   int ID, shapeID, ghs3dShapeID;
 829   int IdShapeRef = 1;
 830   int compoundID =
 831     nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
 832
 833   int *tab, *tabID, *nodeID, *nodeAssigne;
 834   double *coord;
 835   const SMDS_MeshNode **node;
 836
 837   tab    = new int[3];
 838   //tabID  = new int[nbShape];
 839   nodeID = new int[4];
 840   coord  = new double[3];
 841   node   = new const SMDS_MeshNode*[4];
 842
 843   TopoDS_Shape aSolid;
 844   SMDS_MeshNode * aNewNode;
 845   map <int,const SMDS_MeshNode*>::iterator itOnNode;
 846   SMDS_MeshElement* aTet;
 847 #ifdef _DEBUG_
 848   set<int> shapeIDs;
 849 #endif
 850
 851   // Read the file state
 852   fileStat = fstat(fileOpen, &status);
 853   length   = status.st_size;
 854
 855   // Mapping the result file into memory
 856 #ifdef WNT
 857   HANDLE fd = CreateFile(fileName, GENERIC_READ, FILE_SHARE_READ,
 858                          NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
 859   HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
 860                                         0, (DWORD)length, NULL);
 861   ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
 862 #else
 863   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
 864 #endif
 865   mapPtr = ptr;
 866
 867   ptr      = readMapIntLine(ptr, tab);
 868   tetraPtr = ptr;
 869
 870   nbElems      = tab[0];
 871   nbNodes      = tab[1];
 872   nbInputNodes = tab[2];
 873
 874   nodeAssigne = new int[ nbNodes+1 ];
 875
 876   if (nbShape > 0)
 877     aSolid = tabShape[0];
 878
 879   // Reading the nodeId
 880   for (int i=0; i < 4*nbElems; i++)
 881     nodeId = strtol(ptr, &ptr, 10);
 882
 883   MESSAGE("nbInputNodes: "<<nbInputNodes);
 884   MESSAGE("nbEnforcedVertices: "<<nbEnforcedVertices);
 885   // Reading the nodeCoor and update the nodeMap
 886   for (int iNode=1; iNode <= nbNodes; iNode++) {
 887     for (int iCoor=0; iCoor < 3; iCoor++)
 888       coord[ iCoor ] = strtod(ptr, &ptr);
 889     nodeAssigne[ iNode ] = 1;
 890     if ( iNode > (nbInputNodes-nbEnforcedVertices) ) {
 891       // Creating SMESH nodes
 892       // - for enforced vertices
 893       // - for vertices of forced edges
 894       // - for ghs3d nodes
 895       nodeAssigne[ iNode ] = 0;
 896       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
 897       theGhs3dIdToNodeMap.insert(theGhs3dIdToNodeMap.end(), make_pair( iNode, aNewNode ));
 898     }
 899   }
 900
 901   // Reading the number of triangles which corresponds to the number of sub-domains
 902   nbTriangle = strtol(ptr, &ptr, 10);
 903
 904   tabID = new int[nbTriangle];
 905   for (int i=0; i < nbTriangle; i++) {
 906     tabID[i] = 0;
 907     // find the solid corresponding to GHS3D sub-domain following
 908     // the technique proposed in GHS3D manual in chapter
 909     // "B.4 Subdomain (sub-region) assignment"
 910     int nodeId1 = strtol(ptr, &ptr, 10);
 911     int nodeId2 = strtol(ptr, &ptr, 10);
 912     int nodeId3 = strtol(ptr, &ptr, 10);
 913     if ( nbTriangle > 1 ) {
 914       const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
 915       const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
 916       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
 917       try {
 918         OCC_CATCH_SIGNALS;
 919         tabID[i] = findShapeID( theMesh, n1, n2, n3, toMeshHoles );
 920         // -- 0020330: Pb with ghs3d as a submesh
 921         // check that found shape is to be meshed
 922         if ( tabID[i] > 0 ) {
 923           const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] );
 924           bool isToBeMeshed = false;
 925           for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS )
 926             isToBeMeshed = foundShape.IsSame( tabShape[ iS ]);
 927           if ( !isToBeMeshed )
 928             tabID[i] = HOLE_ID;
 929         }
 930         // END -- 0020330: Pb with ghs3d as a submesh
 931 #ifdef _DEBUG_
 932         std::cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << std::endl;
 933 #endif
 934       }
 935       catch ( Standard_Failure & ex)
 936       {
 937 #ifdef _DEBUG_
 938         std::cout << i+1 << " subdomain: Exception caugt: " << ex.GetMessageString() << std::endl;
 939 #endif
 940       }
 941       catch (...) {
 942 #ifdef _DEBUG_
 943         std::cout << i+1 << " subdomain: unknown exception caught " << std::endl;
 944 #endif
 945       }
 946     }
 947   }
 948
 949   shapePtr = ptr;
 950
 951   if ( nbTriangle <= nbShape ) // no holes
 952     toMeshHoles = true; // not avoid creating tetras in holes
 953
 954   // Associating the tetrahedrons to the shapes
 955   shapeID = compoundID;
 956   for (int iElem = 0; iElem < nbElems; iElem++) {
 957     for (int iNode = 0; iNode < 4; iNode++) {
 958       ID = strtol(tetraPtr, &tetraPtr, 10);
 959       itOnNode = theGhs3dIdToNodeMap.find(ID);
 960       node[ iNode ] = itOnNode->second;
 961       nodeID[ iNode ] = ID;
 962     }
 963     // We always run GHS3D with "to mesh holes"==TRUE but we must not create
 964     // tetras within holes depending on hypo option,
 965     // so we first check if aTet is inside a hole and then create it
 966     //aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
 967     if ( nbTriangle > 1 ) {
 968       shapeID = HOLE_ID; // negative shapeID means not to create tetras if !toMeshHoles
 969       ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef;
 970       if ( tabID[ ghs3dShapeID ] == 0 ) {
 971         TopAbs_State state;
 972         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
 973         if ( toMeshHoles || state == TopAbs_IN )
 974           shapeID = theMeshDS->ShapeToIndex( aSolid );
 975         tabID[ ghs3dShapeID ] = shapeID;
 976       }
 977       else
 978         shapeID = tabID[ ghs3dShapeID ];
 979     }
 980     else if ( nbShape > 1 ) {
 981       // Case where nbTriangle == 1 while nbShape == 2 encountered
 982       // with compound of 2 boxes and "To mesh holes"==False,
 983       // so there are no subdomains specified for each tetrahedron.
 984       // Try to guess a solid by a node already bound to shape
 985       shapeID = 0;
 986       for ( int i=0; i<4 && shapeID==0; i++ ) {
 987         if ( nodeAssigne[ nodeID[i] ] == 1 &&
 988              node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
 989              node[i]->getshapeId() > 1 )
 990         {
 991           shapeID = node[i]->getshapeId();
 992         }
 993       }
 994       if ( shapeID==0 ) {
 995         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
 996         shapeID = theMeshDS->ShapeToIndex( aSolid );
 997       }
 998     }
 999     // set new nodes and tetrahedron onto the shape
1000     for ( int i=0; i<4; i++ ) {
1001       if ( nodeAssigne[ nodeID[i] ] == 0 ) {
1002         if ( shapeID != HOLE_ID )
1003           theMeshDS->SetNodeInVolume( node[i], shapeID );
1004         nodeAssigne[ nodeID[i] ] = shapeID;
1005       }
1006     }
1007     if ( toMeshHoles || shapeID != HOLE_ID ) {
1008       aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
1009       theMeshDS->SetMeshElementOnShape( aTet, shapeID );
1010     }
1011 #ifdef _DEBUG_
1012     shapeIDs.insert( shapeID );
1013 #endif
1014   }
1015   // Remove nodes of tetras inside holes if !toMeshHoles
1016   if ( !toMeshHoles ) {
1017     itOnNode = theGhs3dIdToNodeMap.find( nbInputNodes );
1018     for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
1019       ID = itOnNode->first;
1020       if ( nodeAssigne[ ID ] == HOLE_ID )
1021         theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
1022     }
1023   }
1024
1025   if ( nbElems )
1026     cout << nbElems << " tetrahedrons have been associated to " << nbShape << " shapes" << endl;
1027 #ifdef WNT
1028   UnmapViewOfFile(mapPtr);
1029   CloseHandle(hMapObject);
1030   CloseHandle(fd);
1031 #else
1032   munmap(mapPtr, length);
1033 #endif
1034   close(fileOpen);
1035
1036   delete [] tab;
1037   delete [] tabID;
1038   delete [] nodeID;
1039   delete [] coord;
1040   delete [] node;
1041   delete [] nodeAssigne;
1042
1043 #ifdef _DEBUG_
1044   shapeIDs.erase(-1);
1045   if ( shapeIDs.size() != nbShape ) {
1046     std::cout << "Only " << shapeIDs.size() << " solids of " << nbShape << " found" << std::endl;
1047     for (int i=0; i<nbShape; i++) {
1048       shapeID = theMeshDS->ShapeToIndex( tabShape[i] );
1049       if ( shapeIDs.find( shapeID ) == shapeIDs.end() )
1050         std::cout << "  Solid #" << shapeID << " not found" << std::endl;
1051     }
1052   }
1053 #endif
1054
1055   return true;
1056 }
1057
1058 //=======================================================================
1059 //function : readResultFile
1060 //purpose  :
1061 //=======================================================================
1062
1063 static bool readResultFile(const int                      fileOpen,
1064 #ifdef WNT
1065                            const char*                    fileName,
1066 #endif
1067                            SMESH_Mesh&                    theMesh,
1068                            TopoDS_Shape                   aSolid,
1069                            vector <const SMDS_MeshNode*>& theNodeByGhs3dId,
1070                            int                            nbEnforcedVertices)
1071 {
1072   SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
1073
1074   Kernel_Utils::Localizer loc;
1075   struct stat  status;
1076   size_t       length;
1077
1078   char *ptr, *mapPtr;
1079   char *tetraPtr;
1080   char *shapePtr;
1081
1082   int fileStat;
1083   int nbElems, nbNodes, nbInputNodes;
1084   int nodeId, triangleId;
1085   int nbTriangle;
1086   int ID, shapeID;
1087
1088   int *tab;
1089   double *coord;
1090   const SMDS_MeshNode **node;
1091
1092   tab   = new int[3];
1093   coord = new double[3];
1094   node  = new const SMDS_MeshNode*[4];
1095
1096   SMDS_MeshNode * aNewNode;
1097   map <int,const SMDS_MeshNode*>::iterator IdNode;
1098   SMDS_MeshElement* aTet;
1099
1100   // Read the file state
1101   fileStat = fstat(fileOpen, &status);
1102   length   = status.st_size;
1103
1104   // Mapping the result file into memory
1105 #ifdef WNT
1106   HANDLE fd = CreateFile(fileName, GENERIC_READ, FILE_SHARE_READ,
1107                          NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
1108   HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
1109                                         0, (DWORD)length, NULL);
1110   ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
1111 #else
1112   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
1113 #endif
1114   mapPtr = ptr;
1115
1116   ptr      = readMapIntLine(ptr, tab);
1117   tetraPtr = ptr;
1118
1119   nbElems      = tab[0];
1120   nbNodes      = tab[1];
1121   nbInputNodes = tab[2];
1122
1123   theNodeByGhs3dId.resize( nbNodes );
1124
1125   // Reading the nodeId
1126   for (int i=0; i < 4*nbElems; i++)
1127     nodeId = strtol(ptr, &ptr, 10);
1128
1129   // Issue 0020682. Avoid creating nodes and tetras at place where
1130   // volumic elements already exist
1131   SMESH_ElementSearcher* elemSearcher = 0;
1132   vector< const SMDS_MeshElement* > foundVolumes;
1133   if ( theMesh.NbVolumes() > 0 )
1134     elemSearcher = SMESH_MeshEditor( &theMesh ).GetElementSearcher();
1135
1136   // Reading the nodeCoord and update the nodeMap
1137   shapeID = theMeshDS->ShapeToIndex( aSolid );
1138   for (int iNode=0; iNode < nbNodes; iNode++) {
1139     for (int iCoor=0; iCoor < 3; iCoor++)
1140       coord[ iCoor ] = strtod(ptr, &ptr);
1141     if ((iNode+1) > (nbInputNodes-nbEnforcedVertices)) {
1142       // Issue 0020682. Avoid creating nodes and tetras at place where
1143       // volumic elements already exist
1144       if ( elemSearcher &&
1145            elemSearcher->FindElementsByPoint( gp_Pnt(coord[0],coord[1],coord[2]),
1146                                               SMDSAbs_Volume, foundVolumes ))
1147       {
1148         theNodeByGhs3dId[ iNode ] = 0;
1149       }
1150       else
1151       {
1152         aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
1153         theMeshDS->SetNodeInVolume( aNewNode, shapeID );
1154         theNodeByGhs3dId[ iNode ] = aNewNode;
1155       }
1156     }
1157   }
1158
1159   // Reading the triangles
1160   nbTriangle = strtol(ptr, &ptr, 10);
1161
1162   for (int i=0; i < 3*nbTriangle; i++)
1163     triangleId = strtol(ptr, &ptr, 10);
1164
1165   shapePtr = ptr;
1166
1167   // Associating the tetrahedrons to the shapes
1168   for (int iElem = 0; iElem < nbElems; iElem++) {
1169     for (int iNode = 0; iNode < 4; iNode++) {
1170       ID = strtol(tetraPtr, &tetraPtr, 10);
1171       node[ iNode ] = theNodeByGhs3dId[ ID-1 ];
1172     }
1173     if ( elemSearcher )
1174     {
1175       // Issue 0020682. Avoid creating nodes and tetras at place where
1176       // volumic elements already exist
1177       if ( !node[1] || !node[0] || !node[2] || !node[3] )
1178         continue;
1179       if ( elemSearcher->FindElementsByPoint(( SMESH_MeshEditor::TNodeXYZ(node[0]) +
1180                                                SMESH_MeshEditor::TNodeXYZ(node[1]) +
1181                                                SMESH_MeshEditor::TNodeXYZ(node[2]) +
1182                                                SMESH_MeshEditor::TNodeXYZ(node[3]) ) / 4.,
1183                                              SMDSAbs_Volume, foundVolumes ))
1184         continue;
1185     }
1186     aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
1187     shapeID = theMeshDS->ShapeToIndex( aSolid );
1188     theMeshDS->SetMeshElementOnShape( aTet, shapeID );
1189   }
1190   if ( nbElems )
1191     cout << nbElems << " tetrahedrons have been associated to " << nbTriangle << " shapes" << endl;
1192 #ifdef WNT
1193   UnmapViewOfFile(mapPtr);
1194   CloseHandle(hMapObject);
1195   CloseHandle(fd);
1196 #else
1197   munmap(mapPtr, length);
1198 #endif
1199   close(fileOpen);
1200
1201   delete [] tab;
1202   delete [] coord;
1203   delete [] node;
1204
1205   return true;
1206 }
1207
1208 //=============================================================================
1209 /*!
1210  *Here we are going to use the GHS3D mesher
1211  */
1212 //=============================================================================
1213
1214 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
1215                                 const TopoDS_Shape& theShape)
1216 {
1217   bool Ok(false);
1218   //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
1219
1220   // we count the number of shapes
1221   // _nbShape = countShape( meshDS, TopAbs_SOLID ); -- 0020330: Pb with ghs3d as a submesh
1222   _nbShape = 0;
1223   TopExp_Explorer expBox ( theShape, TopAbs_SOLID );
1224   for ( ; expBox.More(); expBox.Next() )
1225     _nbShape++;
1226
1227   // create bounding box for every shape inside the compound
1228
1229   int iShape = 0;
1230   TopoDS_Shape* tabShape;
1231   double**      tabBox;
1232   tabShape = new TopoDS_Shape[_nbShape];
1233   tabBox   = new double*[_nbShape];
1234   for (int i=0; i<_nbShape; i++)
1235     tabBox[i] = new double[6];
1236   Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
1237
1238   for (expBox.ReInit(); expBox.More(); expBox.Next()) {
1239     tabShape[iShape] = expBox.Current();
1240     Bnd_Box BoundingBox;
1241     BRepBndLib::Add(expBox.Current(), BoundingBox);
1242     BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
1243     tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
1244     tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
1245     tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
1246     iShape++;
1247   }
1248
1249   // a unique working file name
1250   // to avoid access to the same files by eg different users
1251   TCollection_AsciiString aGenericName
1252     = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str();
1253
1254   TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
1255   TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
1256   aFacesFileName  = aGenericName + ".faces";  // in faces
1257   aPointsFileName = aGenericName + ".points"; // in points
1258   aResultFileName = aGenericName + ".noboite";// out points and volumes
1259   aBadResFileName = aGenericName + ".boite";  // out bad result
1260   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
1261   aLogFileName    = aGenericName + ".log";    // log
1262
1263   // -----------------
1264   // make input files
1265   // -----------------
1266
1267   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
1268   ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
1269
1270   Ok =
1271     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
1272   if (!Ok) {
1273     INFOS( "Can't write into " << aFacesFileName);
1274     return error(SMESH_Comment("Can't write into ") << aFacesFileName);
1275   }
1276   map <int,int> aSmdsToGhs3dIdMap;
1277   map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
1278   map<vector<double>,double> enforcedVertices;
1279   int nbEnforcedVertices = 0;
1280   try {
1281     enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
1282     nbEnforcedVertices = enforcedVertices.size();
1283   }
1284   catch(...) {
1285   }
1286
1287   SMESH_MesherHelper helper( theMesh );
1288   helper.SetSubShape( theShape );
1289
1290   {
1291     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
1292
1293     // make prisms on quadrangles
1294     if ( theMesh.NbQuadrangles() > 0 )
1295     {
1296       vector<SMESH_ProxyMesh::Ptr> components;
1297       for (expBox.ReInit(); expBox.More(); expBox.Next())
1298       {
1299         if ( _viscousLayersHyp )
1300         {
1301           proxyMesh = _viscousLayersHyp->Compute( theMesh, expBox.Current() );
1302           if ( !proxyMesh )
1303             return false;
1304         }
1305         StdMeshers_QuadToTriaAdaptor* q2t = new StdMeshers_QuadToTriaAdaptor;
1306         q2t->Compute( theMesh, expBox.Current(), proxyMesh.get() );
1307         components.push_back( SMESH_ProxyMesh::Ptr( q2t ));
1308       }
1309       proxyMesh.reset( new SMESH_ProxyMesh( components ));
1310     }
1311     // build viscous layers
1312     else if ( _viscousLayersHyp )
1313     {
1314       proxyMesh = _viscousLayersHyp->Compute( theMesh, theShape );
1315       if ( !proxyMesh )
1316         return false;
1317     }
1318
1319     Ok = (writePoints( aPointsFile, helper, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, enforcedVertices)
1320           &&
1321           writeFaces ( aFacesFile, *proxyMesh, theShape, aSmdsToGhs3dIdMap ));
1322   }
1323
1324   // Write aSmdsToGhs3dIdMap to temp file
1325   TCollection_AsciiString aSmdsToGhs3dIdMapFileName;
1326   aSmdsToGhs3dIdMapFileName = aGenericName + ".ids";  // ids relation
1327   ofstream aIdsFile  ( aSmdsToGhs3dIdMapFileName.ToCString()  , ios::out);
1328   Ok =
1329     aIdsFile.rdbuf()->is_open();
1330   if (!Ok) {
1331     INFOS( "Can't write into " << aSmdsToGhs3dIdMapFileName);
1332     return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName);
1333   }
1334   aIdsFile << "Smds Ghs3d" << std::endl;
1335   map <int,int>::const_iterator myit;
1336   for (myit=aSmdsToGhs3dIdMap.begin() ; myit != aSmdsToGhs3dIdMap.end() ; ++myit) {
1337     aIdsFile << myit->first << " " << myit->second << std::endl;
1338   }
1339
1340   aFacesFile.close();
1341   aPointsFile.close();
1342   aIdsFile.close();
1343
1344   if ( ! Ok ) {
1345     if ( !_keepFiles ) {
1346       removeFile( aFacesFileName );
1347       removeFile( aPointsFileName );
1348       removeFile( aSmdsToGhs3dIdMapFileName );
1349     }
1350     return error(COMPERR_BAD_INPUT_MESH);
1351   }
1352   removeFile( aResultFileName ); // needed for boundary recovery module usage
1353
1354   // -----------------
1355   // run ghs3d mesher
1356   // -----------------
1357
1358   TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
1359   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
1360   cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
1361
1362   std::cout << std::endl;
1363   std::cout << "Ghs3d execution..." << std::endl;
1364   std::cout << cmd << std::endl;
1365
1366   system( cmd.ToCString() ); // run
1367
1368   std::cout << std::endl;
1369   std::cout << "End of Ghs3d execution !" << std::endl;
1370
1371   // --------------
1372   // read a result
1373   // --------------
1374
1375   // Mapping the result file
1376
1377   int fileOpen;
1378   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
1379   if ( fileOpen < 0 ) {
1380     std::cout << std::endl;
1381     std::cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << std::endl;
1382     std::cout << "Log: " << aLogFileName << std::endl;
1383     Ok = false;
1384   }
1385   else {
1386     bool toMeshHoles =
1387       _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
1388     Ok = readResultFile( fileOpen,
1389 #ifdef WNT
1390                          aResultFileName.ToCString(),
1391 #endif
1392                          theMesh, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
1393                          toMeshHoles, nbEnforcedVertices );
1394   }
1395
1396   // ---------------------
1397   // remove working files
1398   // ---------------------
1399
1400   if ( Ok )
1401   {
1402     if ( !_keepFiles )
1403       removeFile( aLogFileName );
1404   }
1405   else if ( OSD_File( aLogFileName ).Size() > 0 )
1406   {
1407     // get problem description from the log file
1408     _Ghs2smdsConvertor conv( aGhs3dIdToNodeMap );
1409     storeErrorDescription( aLogFileName, conv );
1410   }
1411   else
1412   {
1413     // the log file is empty
1414     removeFile( aLogFileName );
1415     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
1416     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
1417   }
1418
1419   if ( !_keepFiles ) {
1420     removeFile( aFacesFileName );
1421     removeFile( aPointsFileName );
1422     removeFile( aResultFileName );
1423     removeFile( aBadResFileName );
1424     removeFile( aBbResFileName );
1425     removeFile( aSmdsToGhs3dIdMapFileName );
1426   }
1427   std::cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
1428   if ( !Ok )
1429     std::cout << "not ";
1430   std::cout << "treated !" << std::endl;
1431   std::cout << std::endl;
1432
1433   _nbShape = 0;    // re-initializing _nbShape for the next Compute() method call
1434   delete [] tabShape;
1435   delete [] tabBox;
1436
1437   return Ok;
1438 }
1439
1440 //=============================================================================
1441 /*!
1442  *Here we are going to use the GHS3D mesher w/o geometry
1443  */
1444 //=============================================================================
1445 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
1446                                 SMESH_MesherHelper* aHelper)
1447 {
1448   MESSAGE("GHS3DPlugin_GHS3D::Compute()");
1449
1450   //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
1451   TopoDS_Shape theShape = aHelper->GetSubShape();
1452
1453   // a unique working file name
1454   // to avoid access to the same files by eg different users
1455   TCollection_AsciiString aGenericName
1456     = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str();
1457
1458   TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
1459   TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
1460   aFacesFileName  = aGenericName + ".faces";  // in faces
1461   aPointsFileName = aGenericName + ".points"; // in points
1462   aResultFileName = aGenericName + ".noboite";// out points and volumes
1463   aBadResFileName = aGenericName + ".boite";  // out bad result
1464   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
1465   aLogFileName    = aGenericName + ".log";    // log
1466
1467   // -----------------
1468   // make input files
1469   // -----------------
1470
1471   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
1472   ofstream aPointsFile  ( aPointsFileName.ToCString()  , ios::out);
1473   bool Ok =
1474     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
1475
1476   if (!Ok)
1477     return error( SMESH_Comment("Can't write into ") << aPointsFileName);
1478
1479   GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices;
1480   int nbEnforcedVertices = 0;
1481   try {
1482     enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
1483     nbEnforcedVertices = enforcedVertices.size();
1484   }
1485   catch(...) {
1486   }
1487
1488   vector <const SMDS_MeshNode*> aNodeByGhs3dId;
1489
1490   SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
1491   if ( theMesh.NbQuadrangles() > 0 )
1492   {
1493     StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor;
1494     aQuad2Trias->Compute( theMesh );
1495     proxyMesh.reset( aQuad2Trias );
1496   }
1497
1498   Ok = (writeFaces ( aFacesFile, *proxyMesh, aNodeByGhs3dId ) &&
1499         writePoints( aPointsFile, &theMesh, aNodeByGhs3dId,enforcedVertices));
1500
1501   aFacesFile.close();
1502   aPointsFile.close();
1503
1504   if ( ! Ok ) {
1505     if ( !_keepFiles ) {
1506       removeFile( aFacesFileName );
1507       removeFile( aPointsFileName );
1508     }
1509     return error(COMPERR_BAD_INPUT_MESH);
1510   }
1511   removeFile( aResultFileName ); // needed for boundary recovery module usage
1512
1513   // -----------------
1514   // run ghs3d mesher
1515   // -----------------
1516
1517   TCollection_AsciiString cmd =
1518     (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str();
1519   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
1520   cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
1521
1522   system( cmd.ToCString() ); // run
1523
1524   // --------------
1525   // read a result
1526   // --------------
1527   int fileOpen;
1528   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
1529   if ( fileOpen < 0 ) {
1530     std::cout << std::endl;
1531     std::cout << "Error when opening the " << aResultFileName.ToCString() << " file" << std::endl;
1532     std::cout << "Log: " << aLogFileName << std::endl;
1533     std::cout << std::endl;
1534     Ok = false;
1535   }
1536   else {
1537     Ok = readResultFile( fileOpen,
1538 #ifdef WNT
1539                          aResultFileName.ToCString(),
1540 #endif
1541                          theMesh, theShape ,aNodeByGhs3dId, nbEnforcedVertices );
1542   }
1543
1544   // ---------------------
1545   // remove working files
1546   // ---------------------
1547
1548   if ( Ok )
1549   {
1550     if ( !_keepFiles )
1551       removeFile( aLogFileName );
1552   }
1553   else if ( OSD_File( aLogFileName ).Size() > 0 )
1554   {
1555     // get problem description from the log file
1556     _Ghs2smdsConvertor conv( aNodeByGhs3dId );
1557     storeErrorDescription( aLogFileName, conv );
1558   }
1559   else {
1560     // the log file is empty
1561     removeFile( aLogFileName );
1562     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
1563     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
1564   }
1565
1566   if ( !_keepFiles )
1567   {
1568     removeFile( aFacesFileName );
1569     removeFile( aPointsFileName );
1570     removeFile( aResultFileName );
1571     removeFile( aBadResFileName );
1572     removeFile( aBbResFileName );
1573   }
1574
1575   return Ok;
1576 }
1577
1578 //================================================================================
1579 /*!
1580  * \brief Provide human readable text by error code reported by ghs3d
1581  */
1582 //================================================================================
1583
1584 static string translateError(const int errNum)
1585 {
1586   switch ( errNum ) {
1587   case 0:
1588     return "The surface mesh includes a face of type other than edge, "
1589       "triangle or quadrilateral. This face type is not supported.";
1590   case 1:
1591     return "Not enough memory for the face table.";
1592   case 2:
1593     return "Not enough memory.";
1594   case 3:
1595     return "Not enough memory.";
1596   case 4:
1597     return "Face is ignored.";
1598   case 5:
1599     return "End of file. Some data are missing in the file.";
1600   case 6:
1601     return "Read error on the file. There are wrong data in the file.";
1602   case 7:
1603     return "the metric file is inadequate (dimension other than 3).";
1604   case 8:
1605     return "the metric file is inadequate (values not per vertices).";
1606   case 9:
1607     return "the metric file contains more than one field.";
1608   case 10:
1609     return "the number of values in the \".bb\" (metric file) is incompatible with the expected"
1610       "value of number of mesh vertices in the \".noboite\" file.";
1611   case 12:
1612     return "Too many sub-domains.";
1613   case 13:
1614     return "the number of vertices is negative or null.";
1615   case 14:
1616     return "the number of faces is negative or null.";
1617   case 15:
1618     return "A face has a null vertex.";
1619   case 22:
1620     return "incompatible data.";
1621   case 131:
1622     return "the number of vertices is negative or null.";
1623   case 132:
1624     return "the number of vertices is negative or null (in the \".mesh\" file).";
1625   case 133:
1626     return "the number of faces is negative or null.";
1627   case 1000:
1628     return "A face appears more than once in the input surface mesh.";
1629   case 1001:
1630     return "An edge appears more than once in the input surface mesh.";
1631   case 1002:
1632     return "A face has a vertex negative or null.";
1633   case 1003:
1634     return "NOT ENOUGH MEMORY.";
1635   case 2000:
1636     return "Not enough available memory.";
1637   case 2002:
1638     return "Some initial points cannot be inserted. The surface mesh is probably very bad "
1639       "in terms of quality or the input list of points is wrong.";
1640   case 2003:
1641     return "Some vertices are too close to one another or coincident.";
1642   case 2004:
1643     return "Some vertices are too close to one another or coincident.";
1644   case 2012:
1645     return "A vertex cannot be inserted.";
1646   case 2014:
1647     return "There are at least two points considered as coincident.";
1648   case 2103:
1649     return "Some vertices are too close to one another or coincident.";
1650   case 3000:
1651     return "The surface mesh regeneration step has failed.";
1652   case 3009:
1653     return "Constrained edge cannot be enforced.";
1654   case 3019:
1655     return "Constrained face cannot be enforced.";
1656   case 3029:
1657     return "Missing faces.";
1658   case 3100:
1659     return "No guess to start the definition of the connected component(s).";
1660   case 3101:
1661     return "The surface mesh includes at least one hole. The domain is not well defined.";
1662   case 3102:
1663     return "Impossible to define a component.";
1664   case 3103:
1665     return "The surface edge intersects another surface edge.";
1666   case 3104:
1667     return "The surface edge intersects the surface face.";
1668   case 3105:
1669     return "One boundary point lies within a surface face.";
1670   case 3106:
1671     return "One surface edge intersects a surface face.";
1672   case 3107:
1673     return "One boundary point lies within a surface edge.";
1674   case 3108:
1675     return "Insufficient memory ressources detected due to a bad quality surface mesh leading "
1676       "to too many swaps.";
1677   case 3109:
1678     return "Edge is unique (i.e., bounds a hole in the surface).";
1679   case 3122:
1680     return "Presumably, the surface mesh is not compatible with the domain being processed.";
1681   case 3123:
1682     return "Too many components, too many sub-domain.";
1683   case 3209:
1684     return "The surface mesh includes at least one hole. "
1685       "Therefore there is no domain properly defined.";
1686   case 3300:
1687     return "Statistics.";
1688   case 3400:
1689     return "Statistics.";
1690   case 3500:
1691     return "Warning, it is dramatically tedious to enforce the boundary items.";
1692   case 4000:
1693     return "Not enough memory at this time, nevertheless, the program continues. "
1694       "The expected mesh will be correct but not really as large as required.";
1695   case 4002:
1696     return "see above error code, resulting quality may be poor.";
1697   case 4003:
1698     return "Not enough memory at this time, nevertheless, the program continues (warning).";
1699   case 8000:
1700     return "Unknown face type.";
1701   case 8005:
1702   case 8006:
1703     return "End of file. Some data are missing in the file.";
1704   case 9000:
1705     return "A too small volume element is detected.";
1706   case 9001:
1707     return "There exists at least a null or negative volume element.";
1708   case 9002:
1709     return "There exist null or negative volume elements.";
1710   case 9003:
1711     return "A too small volume element is detected. A face is considered being degenerated.";
1712   case 9100:
1713     return "Some element is suspected to be very bad shaped or wrong.";
1714   case 9102:
1715     return "A too bad quality face is detected. This face is considered degenerated.";
1716   case 9112:
1717     return "A too bad quality face is detected. This face is degenerated.";
1718   case 9122:
1719     return "Presumably, the surface mesh is not compatible with the domain being processed.";
1720   case 9999:
1721     return "Abnormal error occured, contact hotline.";
1722   case 23600:
1723     return "Not enough memory for the face table.";
1724   case 23601:
1725     return "The algorithm cannot run further. "
1726       "The surface mesh is probably very bad in terms of quality.";
1727   case 23602:
1728     return "Bad vertex number.";
1729   }
1730   return "";
1731 }
1732
1733 //================================================================================
1734 /*!
1735  * \brief Retrieve from a string given number of integers
1736  */
1737 //================================================================================
1738
1739 static char* getIds( char* ptr, int nbIds, vector<int>& ids )
1740 {
1741   ids.clear();
1742   ids.reserve( nbIds );
1743   while ( nbIds )
1744   {
1745     while ( !isdigit( *ptr )) ++ptr;
1746     if ( ptr[-1] == '-' ) --ptr;
1747     ids.push_back( strtol( ptr, &ptr, 10 ));
1748     --nbIds;
1749   }
1750   return ptr;
1751 }
1752
1753 //================================================================================
1754 /*!
1755  * \brief Retrieve problem description form a log file
1756  *  \retval bool - always false
1757  */
1758 //================================================================================
1759
1760 bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& logFile,
1761                                               const _Ghs2smdsConvertor &     toSmdsConvertor )
1762 {
1763   // open file
1764 #ifdef WNT
1765   int file = ::_open (logFile.ToCString(), _O_RDONLY|_O_BINARY);
1766 #else
1767   int file = ::open (logFile.ToCString(), O_RDONLY);
1768 #endif
1769   if ( file < 0 )
1770     return error( SMESH_Comment("See ") << logFile << " for problem description");
1771
1772   // get file size
1773 //   struct stat status;
1774 //   fstat(file, &status);
1775 //   size_t length = status.st_size;
1776   off_t length = lseek( file, 0, SEEK_END);
1777   lseek( file, 0, SEEK_SET);
1778
1779   // read file
1780   vector< char > buf( length );
1781   int nBytesRead = ::read (file, & buf[0], length);
1782   ::close (file);
1783   char* ptr = & buf[0];
1784   char* bufEnd = ptr + nBytesRead;
1785
1786   SMESH_Comment errDescription;
1787
1788   enum { NODE = 1, EDGE, TRIA, VOL, ID = 1 };
1789
1790   // look for errors "ERR #"
1791
1792   set<string> foundErrorStr; // to avoid reporting same error several times
1793   set<int>    elemErrorNums; // not to report different types of errors with bad elements
1794   while ( ++ptr < bufEnd )
1795   {
1796     if ( strncmp( ptr, "ERR ", 4 ) != 0 )
1797       continue;
1798
1799     list<const SMDS_MeshElement*> badElems;
1800     vector<int> nodeIds;
1801
1802     ptr += 4;
1803     char* errBeg = ptr;
1804     int   errNum = strtol(ptr, &ptr, 10);
1805     switch ( errNum ) { // we treat errors enumerated in [SALOME platform 0019316] issue
1806     case 0015:
1807       // The face number (numfac) with vertices (f 1, f 2, f 3) has a null vertex.
1808       ptr = getIds(ptr, NODE, nodeIds);
1809       ptr = getIds(ptr, TRIA, nodeIds);
1810       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1811       break;
1812     case 1000: // ERR  1000 :  1 3 2
1813       // Face (f 1, f 2, f 3) appears more than once in the input surface mesh.
1814       ptr = getIds(ptr, TRIA, nodeIds);
1815       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1816       break;
1817     case 1001:
1818       // Edge (e1, e2) appears more than once in the input surface mesh
1819       ptr = getIds(ptr, EDGE, nodeIds);
1820       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1821       break;
1822     case 1002:
1823       // Face (f 1, f 2, f 3) has a vertex negative or null
1824       ptr = getIds(ptr, TRIA, nodeIds);
1825       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1826       break;
1827     case 2004:
1828       // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
1829       ptr = getIds(ptr, NODE, nodeIds);
1830       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1831       ptr = getIds(ptr, NODE, nodeIds);
1832       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1833       break;
1834     case 2012:
1835       // Vertex v1 cannot be inserted (warning).
1836       ptr = getIds(ptr, NODE, nodeIds);
1837       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1838       break;
1839     case 2014:
1840       // There are at least two points whose distance is dist, i.e., considered as coincident
1841     case 2103: // ERR  2103 :  16 WITH  3
1842       // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
1843       ptr = getIds(ptr, NODE, nodeIds);
1844       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1845       ptr = getIds(ptr, NODE, nodeIds);
1846       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1847       break;
1848     case 3009:
1849       // Constrained edge (e1, e2) cannot be enforced (warning).
1850       ptr = getIds(ptr, EDGE, nodeIds);
1851       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1852       break;
1853     case 3019:
1854       // Constrained face (f 1, f 2, f 3) cannot be enforced
1855       ptr = getIds(ptr, TRIA, nodeIds);
1856       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1857       break;
1858     case 3103: // ERR  3103 :  1 2 WITH  7 3
1859       // The surface edge (e1, e2) intersects another surface edge (e3, e4)
1860       ptr = getIds(ptr, EDGE, nodeIds);
1861       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1862       ptr = getIds(ptr, EDGE, nodeIds);
1863       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1864       break;
1865     case 3104: // ERR  3104 :  9 10 WITH  1 2 3
1866       // The surface edge (e1, e2) intersects the surface face (f 1, f 2, f 3)
1867       ptr = getIds(ptr, EDGE, nodeIds);
1868       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1869       ptr = getIds(ptr, TRIA, nodeIds);
1870       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1871       break;
1872     case 3105: // ERR  3105 :  8 IN  2 3 5
1873       // One boundary point (say p1) lies within a surface face (f 1, f 2, f 3)
1874       ptr = getIds(ptr, NODE, nodeIds);
1875       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1876       ptr = getIds(ptr, TRIA, nodeIds);
1877       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1878       break;
1879     case 3106:
1880       // One surface edge (say e1, e2) intersects a surface face (f 1, f 2, f 3)
1881       ptr = getIds(ptr, EDGE, nodeIds);
1882       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1883       ptr = getIds(ptr, TRIA, nodeIds);
1884       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1885       break;
1886     case 3107: // ERR  3107 :  2 IN  4 1
1887       // One boundary point (say p1) lies within a surface edge (e1, e2) (stop).
1888       ptr = getIds(ptr, NODE, nodeIds);
1889       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1890       ptr = getIds(ptr, EDGE, nodeIds);
1891       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1892       break;
1893     case 3109: // ERR  3109 :  EDGE  5 6 UNIQUE
1894       // Edge (e1, e2) is unique (i.e., bounds a hole in the surface)
1895       ptr = getIds(ptr, EDGE, nodeIds);
1896       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1897       break;
1898     case 9000: // ERR  9000
1899       //  ELEMENT  261 WITH VERTICES :  7 396 -8 242
1900       //  VOLUME   : -1.11325045E+11 W.R.T. EPSILON   0.
1901       // A too small volume element is detected. Are reported the index of the element,
1902       // its four vertex indices, its volume and the tolerance threshold value
1903       ptr = getIds(ptr, ID, nodeIds);
1904       ptr = getIds(ptr, VOL, nodeIds);
1905       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1906       // even if all nodes found, volume it most probably invisible,
1907       // add its faces to demenstrate it anyhow
1908       {
1909         vector<int> faceNodes( nodeIds.begin(), --nodeIds.end() ); // 012
1910         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1911         faceNodes[2] = nodeIds[3]; // 013
1912         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1913         faceNodes[1] = nodeIds[2]; // 023
1914         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1915         faceNodes[0] = nodeIds[1]; // 123
1916         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1917       }
1918       break;
1919     case 9001: // ERR  9001
1920       //  %% NUMBER OF NEGATIVE VOLUME TETS  :  1
1921       //  %% THE LARGEST NEGATIVE TET        :   1.75376581E+11
1922       //  %%  NUMBER OF NULL VOLUME TETS     :  0
1923       // There exists at least a null or negative volume element
1924       break;
1925     case 9002:
1926       // There exist n null or negative volume elements
1927       break;
1928     case 9003:
1929       // A too small volume element is detected
1930       break;
1931     case 9102:
1932       // A too bad quality face is detected. This face is considered degenerated,
1933       // its index, its three vertex indices together with its quality value are reported
1934       break; // same as next
1935     case 9112: // ERR  9112
1936       //  FACE   2 WITH VERTICES :  4 2 5
1937       //  SMALL INRADIUS :   0.
1938       // A too bad quality face is detected. This face is degenerated,
1939       // its index, its three vertex indices together with its inradius are reported
1940       ptr = getIds(ptr, ID, nodeIds);
1941       ptr = getIds(ptr, TRIA, nodeIds);
1942       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1943       // add triangle edges as it most probably has zero area and hence invisible
1944       {
1945         vector<int> edgeNodes(2);
1946         edgeNodes[0] = nodeIds[0]; edgeNodes[1] = nodeIds[1]; // 0-1
1947         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
1948         edgeNodes[1] = nodeIds[2]; // 0-2
1949         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
1950         edgeNodes[0] = nodeIds[1]; // 1-2
1951         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
1952       }
1953       break;
1954     }
1955
1956     bool isNewError = foundErrorStr.insert( string( errBeg, ptr )).second;
1957     if ( !isNewError )
1958       continue; // not to report same error several times
1959
1960 //     const SMDS_MeshElement* nullElem = 0;
1961 //     bool allElemsOk = ( find( badElems.begin(), badElems.end(), nullElem) == badElems.end());
1962
1963 //     if ( allElemsOk && !badElems.empty() && !elemErrorNums.empty() ) {
1964 //       bool oneMoreErrorType = elemErrorNums.insert( errNum ).second;
1965 //       if ( oneMoreErrorType )
1966 //         continue; // not to report different types of errors with bad elements
1967 //     }
1968
1969     // store bad elements
1970     //if ( allElemsOk ) {
1971       list<const SMDS_MeshElement*>::iterator elem = badElems.begin();
1972       for ( ; elem != badElems.end(); ++elem )
1973         addBadInputElement( *elem );
1974       //}
1975
1976     // make error text
1977     string text = translateError( errNum );
1978     if ( errDescription.find( text ) == text.npos ) {
1979       if ( !errDescription.empty() )
1980         errDescription << "\n";
1981       errDescription << text;
1982     }
1983
1984   } // end while
1985
1986   if ( errDescription.empty() ) { // no errors found
1987     char msg[] = "connection to server failed";
1988     if ( search( &buf[0], bufEnd, msg, msg + strlen(msg)) != bufEnd )
1989       errDescription << "Licence problems.";
1990   }
1991
1992   if ( errDescription.empty() )
1993     errDescription << "See " << logFile << " for problem description";
1994   else
1995     errDescription << "\nSee " << logFile << " for more information";
1996
1997   return error( errDescription );
1998 }
1999
2000 //================================================================================
2001 /*!
2002  * \brief Creates _Ghs2smdsConvertor
2003  */
2004 //================================================================================
2005
2006 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const map <int,const SMDS_MeshNode*> & ghs2NodeMap)
2007   :_ghs2NodeMap( & ghs2NodeMap ), _nodeByGhsId( 0 )
2008 {
2009 }
2010
2011 //================================================================================
2012 /*!
2013  * \brief Creates _Ghs2smdsConvertor
2014  */
2015 //================================================================================
2016
2017 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const vector <const SMDS_MeshNode*> &  nodeByGhsId)
2018   : _ghs2NodeMap( 0 ), _nodeByGhsId( &nodeByGhsId )
2019 {
2020 }
2021
2022 //================================================================================
2023 /*!
2024  * \brief Return SMDS element by ids of GHS3D nodes
2025  */
2026 //================================================================================
2027
2028 const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector<int>& ghsNodes) const
2029 {
2030   size_t nbNodes = ghsNodes.size();
2031   vector<const SMDS_MeshNode*> nodes( nbNodes, 0 );
2032   for ( size_t i = 0; i < nbNodes; ++i ) {
2033     int ghsNode = ghsNodes[ i ];
2034     if ( _ghs2NodeMap ) {
2035       map <int,const SMDS_MeshNode*>::const_iterator in = _ghs2NodeMap->find( ghsNode);
2036       if ( in == _ghs2NodeMap->end() )
2037         return 0;
2038       nodes[ i ] = in->second;
2039     }
2040     else {
2041       if ( ghsNode < 1 || ghsNode > _nodeByGhsId->size() )
2042         return 0;
2043       nodes[ i ] = (*_nodeByGhsId)[ ghsNode-1 ];
2044     }
2045   }
2046   if ( nbNodes == 1 )
2047     return nodes[0];
2048
2049   if ( nbNodes == 2 ) {
2050     const SMDS_MeshElement* edge= SMDS_Mesh::FindEdge( nodes[0], nodes[1] );
2051     if ( !edge )
2052       edge = new SMDS_LinearEdge( nodes[0], nodes[1] );
2053     return edge;
2054   }
2055   if ( nbNodes == 3 ) {
2056     const SMDS_MeshElement* face = SMDS_Mesh::FindFace( nodes );
2057     if ( !face )
2058       face = new SMDS_FaceOfNodes( nodes[0], nodes[1], nodes[2] );
2059     return face;
2060   }
2061   if ( nbNodes == 4 )
2062     return new SMDS_VolumeOfNodes( nodes[0], nodes[1], nodes[2], nodes[3] );
2063
2064   return 0;
2065 }
2066
2067
2068 //=============================================================================
2069 /*!
2070  *
2071  */
2072 //=============================================================================
2073 bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh,
2074                                  const TopoDS_Shape& aShape,
2075                                  MapShapeNbElems& aResMap)
2076 {
2077   int nbtri = 0, nbqua = 0;
2078   double fullArea = 0.0;
2079   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
2080     TopoDS_Face F = TopoDS::Face( exp.Current() );
2081     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
2082     MapShapeNbElemsItr anIt = aResMap.find(sm);
2083     if( anIt==aResMap.end() ) {
2084       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2085       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
2086                                             "Submesh can not be evaluated",this));
2087       return false;
2088     }
2089     std::vector<int> aVec = (*anIt).second;
2090     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
2091     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
2092     GProp_GProps G;
2093     BRepGProp::SurfaceProperties(F,G);
2094     double anArea = G.Mass();
2095     fullArea += anArea;
2096   }
2097
2098   // collect info from edges
2099   int nb0d_e = 0, nb1d_e = 0;
2100   bool IsQuadratic = false;
2101   bool IsFirst = true;
2102   TopTools_MapOfShape tmpMap;
2103   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
2104     TopoDS_Edge E = TopoDS::Edge(exp.Current());
2105     if( tmpMap.Contains(E) )
2106       continue;
2107     tmpMap.Add(E);
2108     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
2109     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
2110     std::vector<int> aVec = (*anIt).second;
2111     nb0d_e += aVec[SMDSEntity_Node];
2112     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
2113     if(IsFirst) {
2114       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
2115       IsFirst = false;
2116     }
2117   }
2118   tmpMap.Clear();
2119
2120   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
2121
2122   GProp_GProps G;
2123   BRepGProp::VolumeProperties(aShape,G);
2124   double aVolume = G.Mass();
2125   double tetrVol = 0.1179*ELen*ELen*ELen;
2126   double CoeffQuality = 0.9;
2127   int nbVols = int(aVolume/tetrVol/CoeffQuality);
2128   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
2129   int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
2130   std::vector<int> aVec(SMDSEntity_Last);
2131   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
2132   if( IsQuadratic ) {
2133     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
2134     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
2135     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
2136   }
2137   else {
2138     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
2139     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
2140     aVec[SMDSEntity_Pyramid] = nbqua;
2141   }
2142   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
2143   aResMap.insert(std::make_pair(sm,aVec));
2144
2145   return true;
2146 }
2147