]> SALOME platform Git repositories - plugins/ghs3dplugin.git/commitdiff
Salome HOME
Merge from BR_ENFORCED_MESH
authorgdd <gdd>
Mon, 2 May 2011 16:52:27 +0000 (16:52 +0000)
committergdd <gdd>
Mon, 2 May 2011 16:52:27 +0000 (16:52 +0000)
src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx
src/GHS3DPlugin/GHS3DPlugin_GHS3D.hxx
src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.cxx
src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.hxx
src/GHS3DPlugin/GHS3DPlugin_Hypothesis.cxx
src/GHS3DPlugin/GHS3DPlugin_Hypothesis.hxx
src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.cxx
src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.hxx
src/GHS3DPlugin/Makefile.am
src/GHS3DPlugin/libmesh5.c [new file with mode: 0644]
src/GHS3DPlugin/libmesh5.h [new file with mode: 0755]

index 48caf143285bee60c14820de211edd2c79abc88e..c72664391c9f1343b6c04da6e80c60d8229fc278 100644 (file)
 #include "GHS3DPlugin_GHS3D.hxx"
 #include "GHS3DPlugin_Hypothesis.hxx"
 
-
 #include <Basics_Utils.hxx>
 
-#include "SMESH_Gen.hxx"
+//#include "SMESH_Gen.hxx"
+#include <SMESH_Client.hxx>
 #include "SMESH_Mesh.hxx"
 #include "SMESH_Comment.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_MeshEditor.hxx"
+#include "SMESH_OctreeNode.hxx"
 
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
@@ -46,7 +47,9 @@
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepBndLib.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
 #include <BRepClass3d_SolidClassifier.hxx>
+#include <BRepExtrema_DistShapeShape.hxx>
 #include <BRepTools.hxx>
 #include <BRep_Tool.hxx>
 #include <Bnd_Box.hxx>
@@ -74,6 +77,7 @@
 #else
 #include <sys/sysinfo.h>
 #endif
+#include <algorithm>
 
 using namespace std;
 
@@ -101,6 +105,10 @@ extern "C"
 
 #define HOLE_ID -1
 
+#ifndef GHS3D_VERSION
+#define GHS3D_VERSION 41
+#endif
+
 typedef const list<const SMDS_MeshFace*> TTriaList;
 
 static void removeFile( const TCollection_AsciiString& fileName )
@@ -179,6 +187,7 @@ bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&         aMesh,
   return true;
 }
 
+
 //=======================================================================
 //function : findShape
 //purpose  : 
@@ -238,394 +247,6 @@ static char* readMapIntLine(char* ptr, int tab[]) {
   return ptr;
 }
 
-//=======================================================================
-//function : countShape
-//purpose  :
-//=======================================================================
-
-// template < class Mesh, class Shape >
-// static int countShape( Mesh* mesh, Shape shape ) {
-//   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
-//   int nbShape = 0;
-//   for ( ; expShape.More(); expShape.Next() ) {
-//       nbShape++;
-//   }
-//   return nbShape;
-// }
-
-//=======================================================================
-//function : writeFaces
-//purpose  : 
-//=======================================================================
-
-static bool writeFaces (ofstream &             theFile,
-                        const SMESH_ProxyMesh& theMesh,
-                        const TopoDS_Shape&    theShape,
-                        const map <int,int> &  theSmdsToGhs3dIdMap)
-{
-  // record structure:
-  //
-  // NB_ELEMS DUMMY_INT
-  // Loop from 1 to NB_ELEMS
-  // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
-
-  TopoDS_Shape aShape;
-  const SMESHDS_SubMesh* theSubMesh;
-  const SMDS_MeshElement* aFace;
-  const char* space    = "  ";
-  const int   dummyint = 0;
-  map<int,int>::const_iterator itOnMap;
-  SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
-  int nbNodes, aSmdsID;
-
-  // count triangles bound to geometry
-  int nbTriangles = 0;
-
-  TopTools_IndexedMapOfShape facesMap;
-  TopExp::MapShapes( theShape, TopAbs_FACE, facesMap );
-
-  for ( int i = 1; i <= facesMap.Extent(); ++i )
-    if (( theSubMesh  = theMesh.GetSubMesh( facesMap(i))))
-      nbTriangles += theSubMesh->NbElements();
-
-  std::cout << "    " << facesMap.Extent() << " shapes of 2D dimension" << std::endl;
-  std::cout << std::endl;
-
-  theFile << space << nbTriangles << space << dummyint << std::endl;
-
-  for ( int i = 1; i <= facesMap.Extent(); i++ )
-  {
-    aShape = facesMap(i);
-    theSubMesh = theMesh.GetSubMesh(aShape);
-    if ( !theSubMesh ) continue;
-    itOnSubMesh = theSubMesh->GetElements();
-    while ( itOnSubMesh->more() )
-    {
-      aFace = itOnSubMesh->next();
-      nbNodes = aFace->NbCornerNodes();
-
-      theFile << space << nbNodes;
-
-      itOnSubFace = aFace->nodesIterator();
-      while ( itOnSubFace->more() && nbNodes-- ) {
-        // find GHS3D ID
-        aSmdsID = itOnSubFace->next()->GetID();
-        itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
-        // if ( itOnMap == theSmdsToGhs3dIdMap.end() ) {
-        //   cout << "not found node: " << aSmdsID << endl;
-        //   return false;
-        // }
-        ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
-
-        theFile << space << (*itOnMap).second;
-      }
-
-      // (NB_NODES + 1) times: DUMMY_INT
-      nbNodes = aFace->NbCornerNodes();
-      for ( int j=0; j<=nbNodes; j++)
-        theFile << space << dummyint;
-
-      theFile << std::endl;
-    }
-  }
-
-  return true;
-}
-
-//=======================================================================
-//function : writeFaces
-//purpose  : Write Faces in case if generate 3D mesh w/o geometry
-//=======================================================================
-
-static bool writeFaces (ofstream &                      theFile,
-                        const SMESH_ProxyMesh&          theMesh,
-                        vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
-{
-  // record structure:
-  //
-  // NB_ELEMS DUMMY_INT
-  // Loop from 1 to NB_ELEMS
-  //   NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
-
-  int nbNodes, nbTriangles = 0;
-  map< const SMDS_MeshNode*,int >::iterator it;
-  SMDS_ElemIteratorPtr nodeIt;
-  const SMDS_MeshElement* elem;
-
-  const char* space    = "  ";
-  const int   dummyint = 0;
-
-  // count faces
-  nbTriangles = theMesh.NbFaces();
-
-  if ( nbTriangles == 0 )
-    return false;
-
-  std::cout << "The initial 2D mesh contains " << nbTriangles << " faces and ";
-
-  // NB_ELEMS DUMMY_INT
-  theFile << space << nbTriangles << space << dummyint << std::endl;
-
-
-  map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
-
-  // Loop from 1 to NB_ELEMS
-
-  SMDS_ElemIteratorPtr eIt = theMesh.GetFaces();
-  while ( eIt->more() )
-  {
-    elem = eIt->next();
-    // NB_NODES PER FACE
-    nbNodes = elem->NbCornerNodes();
-    theFile << space << nbNodes;
-
-    // NODE_NB_1 NODE_NB_2 ...
-    nodeIt = elem->nodesIterator();
-    while ( nodeIt->more() && nbNodes-- )
-    {
-      // find GHS3D ID
-      const SMDS_MeshNode* node = castToNode( nodeIt->next() );
-      int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
-      it = aNodeToGhs3dIdMap.insert( make_pair( node, newId )).first;
-      theFile << space << it->second;
-    }
-
-    // (NB_NODES + 1) times: DUMMY_INT
-    nbNodes = elem->NbCornerNodes();
-    for ( int i=0; i<=nbNodes; i++)
-      theFile << space << dummyint;
-    theFile << std::endl;
-  }
-
-  // put nodes to theNodeByGhs3dId vector
-  theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
-  map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
-  for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
-  {
-    theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
-  }
-
-  return true;
-}
-
-//=======================================================================
-//function : writePoints
-//purpose  : 
-//=======================================================================
-
-static bool writePoints (ofstream &                       theFile,
-                         SMESH_MesherHelper&              theHelper,
-                         map <int,int> &                  theSmdsToGhs3dIdMap,
-                         map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap,
-                         map<vector<double>,double> &     theEnforcedVertices)
-{
-  // record structure:
-  //
-  // NB_NODES
-  // Loop from 1 to NB_NODES
-  //   X Y Z DUMMY_INT
-
-  SMESHDS_Mesh * theMesh = theHelper.GetMeshDS();
-  int nbNodes = theMesh->NbNodes();
-  if ( nbNodes == 0 )
-    return false;
-
-  int nbEnforcedVertices = theEnforcedVertices.size();
-
-  int aGhs3dID = 1;
-  SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
-  const SMDS_MeshNode* node;
-
-  // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
-  // The problem is in nodes on degenerated edges, we need to skip them
-  if ( theHelper.HasDegeneratedEdges() )
-  {
-    // here we decrease total nb of nodes by nb of nodes on degenerated edges
-    set<int> checkedSM;
-    for (TopExp_Explorer e(theMesh->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
-    {
-      SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
-      if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() )) {
-        if ( sm->GetSubMeshDS() )
-          nbNodes -= sm->GetSubMeshDS()->NbNodes();
-      }
-    }
-  }
-
-  const bool isQuadMesh = 
-    theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
-    theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
-    theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC );
-  if ( isQuadMesh )
-  {
-    // descrease nbNodes by nb of medium nodes
-    while ( it->more() )
-    {
-      node = it->next();
-      if ( !theHelper.IsDegenShape( node->getshapeId() ))
-        nbNodes -= int( theHelper.IsMedium( node ));
-    }
-    it = theMesh->nodesIterator();
-  }
-
-  const char* space    = "  ";
-  const int   dummyint = 0;
-
-  // NB_NODES
-  std::cout << std::endl;
-  std::cout << "The initial 2D mesh contains :" << std::endl;
-  std::cout << "    " << nbNodes << " nodes" << std::endl;
-  if (nbEnforcedVertices > 0) {
-    std::cout << "    " << nbEnforcedVertices << " enforced vertices" << std::endl;
-  }
-  std::cout << std::endl;
-  std::cout << "Start writing in 'points' file ..." << std::endl;
-  theFile << space << nbNodes << std::endl;
-
-  // Loop from 1 to NB_NODES
-
-  while ( it->more() )
-  {
-    node = it->next();
-    if (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238
-        theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674
-      continue;
-
-    theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
-    theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
-    aGhs3dID++;
-
-    // X Y Z DUMMY_INT
-    theFile
-    << space << node->X()
-    << space << node->Y()
-    << space << node->Z()
-    << space << dummyint;
-
-    theFile << std::endl;
-
-  }
-  
-  // Iterate over the enforced vertices
-  GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
-  const TopoDS_Shape shapeToMesh = theMesh->ShapeToMesh();
-  for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
-    double x = vertexIt->first[0];
-    double y = vertexIt->first[1];
-    double z = vertexIt->first[2];
-    // Test if point is inside shape to mesh
-    gp_Pnt myPoint(x,y,z);
-    BRepClass3d_SolidClassifier scl(shapeToMesh);
-    scl.Perform(myPoint, 1e-7);
-    TopAbs_State result = scl.State();
-    if ( result == TopAbs_IN ) {
-        MESSAGE("Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second);
-        // X Y Z PHY_SIZE DUMMY_INT
-        theFile
-        << space << x
-        << space << y
-        << space << z
-        << space << vertexIt->second
-        << space << dummyint;
-    
-        theFile << std::endl;
-    }
-    else
-        MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added ");
-  }
-  
-  
-  std::cout << std::endl;
-  std::cout << "End writing in 'points' file." << std::endl;
-
-  return true;
-}
-
-//=======================================================================
-//function : writePoints
-//purpose  : 
-//=======================================================================
-
-static bool writePoints (ofstream &                            theFile,
-                         SMESH_Mesh *                          theMesh,
-                         const vector <const SMDS_MeshNode*> & theNodeByGhs3dId,
-                         map<vector<double>,double> &          theEnforcedVertices)
-{
-  // record structure:
-  //
-  // NB_NODES
-  // Loop from 1 to NB_NODES
-  //   X Y Z DUMMY_INT
-  
-  int nbNodes = theNodeByGhs3dId.size();
-  if ( nbNodes == 0 )
-    return false;
-
-  int nbEnforcedVertices = theEnforcedVertices.size();
-
-  const char* space    = "  ";
-  const int   dummyint = 0;
-
-  const SMDS_MeshNode* node;
-
-  // NB_NODES
-  std::cout << std::endl;
-  std::cout << "The initial 2D mesh contains :" << std::endl;
-  std::cout << "    " << nbNodes << " nodes" << std::endl;
-  std::cout << "    " << nbEnforcedVertices << " enforced vertices" << std::endl;
-  std::cout << std::endl;
-  std::cout << "Start writing in 'points' file ..." << std::endl;
-  theFile << space << nbNodes << std::endl;
-
-  // Loop from 1 to NB_NODES
-
-  vector<const SMDS_MeshNode*>::const_iterator nodeIt = theNodeByGhs3dId.begin();
-  vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
-  for ( ; nodeIt != after; ++nodeIt )
-  {
-    node = *nodeIt;
-
-    // X Y Z DUMMY_INT
-    theFile
-    << space << node->X()
-    << space << node->Y()
-    << space << node->Z()
-    << space << dummyint;
-
-    theFile << std::endl;
-
-  }
-  
-  // Iterate over the enforced vertices
-  GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
-  auto_ptr< SMESH_ElementSearcher > pntCls ( SMESH_MeshEditor( theMesh ).GetElementSearcher());
-  for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
-    double x = vertexIt->first[0];
-    double y = vertexIt->first[1];
-    double z = vertexIt->first[2];
-    // Test if point is inside shape to mesh
-    gp_Pnt myPoint(x,y,z);
-    TopAbs_State result = pntCls->GetPointState( myPoint );
-    if ( result == TopAbs_IN ) {
-        std::cout << "Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second << std::endl;
-
-        // X Y Z PHY_SIZE DUMMY_INT
-        theFile
-        << space << x
-        << space << y
-        << space << z
-        << space << vertexIt->second
-        << space << dummyint;
-    
-        theFile << std::endl;
-    }
-  }
-  std::cout << std::endl;
-  std::cout << "End writing in 'points' file." << std::endl;
-
-  return true;
-}
-
 //================================================================================
 /*!
  * \brief returns true if a triangle defined by the nodes is a temporary face on a
@@ -789,34 +410,2029 @@ static int findShapeID(SMESH_Mesh&          mesh,
         case TopAbs_VERTEX: tol = BRep_Tool::Tolerance( TopoDS::Vertex( nodeShape )); break;
         default:;
         }
-      gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
-      BRepAdaptor_Surface surface( geomFace );
-      uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
-      if ( uvOK ) {
-        // normale to geomFace at UV
-        gp_Vec du, dv;
-        surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
-        geomNormal = du ^ dv;
-        if ( geomFace.Orientation() == TopAbs_REVERSED )
-          geomNormal.Reverse();
-        geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
+      gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
+      BRepAdaptor_Surface surface( geomFace );
+      uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
+      if ( uvOK ) {
+        // normale to geomFace at UV
+        gp_Vec du, dv;
+        surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
+        geomNormal = du ^ dv;
+        if ( geomFace.Orientation() == TopAbs_REVERSED )
+          geomNormal.Reverse();
+        geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
+      }
+    }
+  }
+  if ( !geomNormalOK)
+    return invalidID;
+
+  // compare normals
+  bool isReverse = ( meshNormal * geomNormal ) < 0;
+  if ( !isReverse )
+    return meshDS->ShapeToIndex( solid1 );
+
+  if ( solids.Extent() == 1 )
+    return HOLE_ID; // we are inside a hole
+  else
+    return meshDS->ShapeToIndex( solids(2) );
+}
+
+//=======================================================================
+//function : countShape
+//purpose  :
+//=======================================================================
+
+// template < class Mesh, class Shape >
+// static int countShape( Mesh* mesh, Shape shape ) {
+//   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
+//   TopTools_MapOfShape mapShape;
+//   int nbShape = 0;
+//   for ( ; expShape.More(); expShape.Next() ) {
+//     if (mapShape.Add(expShape.Current())) {
+//       nbShape++;
+//     }
+//   }
+//   return nbShape;
+// }
+
+//=======================================================================
+//function : getShape
+//purpose  :
+//=======================================================================
+
+// template < class Mesh, class Shape, class Tab >
+// void getShape(Mesh* mesh, Shape shape, Tab *t_Shape) {
+//   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
+//   TopTools_MapOfShape mapShape;
+//   for ( int i=0; expShape.More(); expShape.Next() ) {
+//     if (mapShape.Add(expShape.Current())) {
+//       t_Shape[i] = expShape.Current();
+//       i++;
+//     }
+//   }
+//   return;
+// }
+
+// //=======================================================================
+// //function : findEdgeID
+// //purpose  :
+// //=======================================================================
+// 
+// static int findEdgeID(const SMDS_MeshNode* aNode,
+//                       const SMESHDS_Mesh*  theMesh,
+//                       const int            nEdge,
+//                       const TopoDS_Shape*  t_Edge) {
+// 
+//   TopoDS_Shape aPntShape, foundEdge;
+//   TopoDS_Vertex aVertex;
+//   gp_Pnt aPnt( aNode->X(), aNode->Y(), aNode->Z() );
+// 
+//   int foundInd, ind;
+//   double nearest = RealLast(), *t_Dist;
+//   double epsilon = Precision::Confusion();
+// 
+//   t_Dist = new double[ nEdge ];
+//   aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape();
+//   aVertex   = TopoDS::Vertex( aPntShape );
+// 
+//   for ( ind=0; ind < nEdge; ind++ ) {
+//     BRepExtrema_DistShapeShape aDistance ( aVertex, t_Edge[ind] );
+//     t_Dist[ind] = aDistance.Value();
+//     if ( t_Dist[ind] < nearest ) {
+//       nearest   = t_Dist[ind];
+//       foundEdge = t_Edge[ind];
+//       foundInd  = ind;
+//       if ( nearest < epsilon )
+//         ind = nEdge;
+//     }
+//   }
+// 
+//   delete [] t_Dist;
+//   return theMesh->ShapeToIndex( foundEdge );
+// }
+
+
+//=======================================================================
+//function : readGMFFile
+//purpose  : read GMF file with geometry associated to mesh
+// TODO
+//=======================================================================
+
+// static bool readGMFFile(
+//                         const int                       fileOpen,
+//                         const char*                     theFileName, 
+//                         SMESH_Mesh&                     theMesh,
+//                         const int                       nbShape,
+//                         const TopoDS_Shape*             tabShape,
+//                         double**                        tabBox,
+//                         map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
+//                         bool                            toMeshHoles,
+//                         int                             nbEnforcedVertices,
+//                         int                             nbEnforcedNodes,
+//                         TIDSortedNodeSet &              theEnforcedNodes,
+//                         TIDSortedElemSet &              theEnforcedTriangles,
+//                         TIDSortedElemSet &              theEnforcedQuadrangles)
+// {
+//   TopoDS_Shape aShape;
+//   TopoDS_Vertex aVertex;
+//   SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
+//   int nbElem = 0, nbRef = 0, IdShapeRef = 1;
+//   int *tabID;
+//   int aGMFNodeID = 0;
+//   int compoundID =
+//     nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
+//   int tetraShapeID = compoundID;
+//   double epsilon = Precision::Confusion();
+//   int *nodeAssigne, *GMFNodeAssigne;
+//   SMDS_MeshNode** GMFNode;
+//   TopoDS_Shape *tabCorner, *tabEdge;
+//   std::map <GmfKwdCod,int> tabRef;
+//   
+//   
+//   int ver, dim;
+//   MESSAGE("Read " << theFileName << " file");
+//   int InpMsh = GmfOpenMesh(theFileName, GmfRead, &ver, &dim);
+//   if (!InpMsh)
+//     return false;
+//   
+//   // ===========================
+//   // Fill the tabID array: BEGIN
+//   // ===========================
+//   
+//   /*
+//   The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
+//   */
+//   Kernel_Utils::Localizer loc;
+//   struct stat status;
+//   size_t      length;
+// 
+//   char *ptr, *mapPtr;
+//   char *tetraPtr;
+//   int *tab = new int[3];
+//   
+//   // Read the file state
+//   fstat(fileOpen, &status);
+//   length   = status.st_size;
+//   
+//   // Mapping the result file into memory
+// #ifdef WNT
+//   HANDLE fd = CreateFile(theFileName, GENERIC_READ, FILE_SHARE_READ,
+//                          NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
+//   HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
+//                                         0, (DWORD)length, NULL);
+//   ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
+// #else
+//   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
+// #endif
+//   mapPtr = ptr;
+// 
+//   ptr      = readMapIntLine(ptr, tab);
+//   tetraPtr = ptr;
+// 
+//   nbElem            = tab[0];
+//   int nbNodes       = tab[1];
+//   
+//   for (int i=0; i < 4*nbElem; i++)
+//     strtol(ptr, &ptr, 10);
+//   
+//   for (int iNode=1; iNode <= nbNodes; iNode++)
+//     for (int iCoor=0; iCoor < 3; iCoor++)
+//       strtod(ptr, &ptr);
+// 
+//     
+//   // Reading the number of triangles which corresponds to the number of sub-domains
+//   int nbTriangle = strtol(ptr, &ptr, 10);
+// 
+//   
+//   // The keyword does not exist yet => to update when it is created
+// //   int nbTriangle = GmfStatKwd(InpMsh, GmfSubdomain);
+// //   int id_tri[3];
+// 
+// 
+//   tabID = new int[nbTriangle];
+//   for (int i=0; i < nbTriangle; i++) {
+//     tabID[i] = 0;
+//     int nodeId1, nodeId2, nodeId3;
+//     // find the solid corresponding to GHS3D sub-domain following
+//     // the technique proposed in GHS3D manual in chapter
+//     // "B.4 Subdomain (sub-region) assignment"
+// 
+//     nodeId1 = strtol(ptr, &ptr, 10);
+//     nodeId2 = strtol(ptr, &ptr, 10);
+//     nodeId3 = strtol(ptr, &ptr, 10);
+// 
+// //   // The keyword does not exist yet => to update when it is created
+// //     GmfGetLin(InpMsh, GmfSubdomain, &id_tri[0], &id_tri[1], &id_tri[2]);
+// //     nodeId1 = id_tri[0];
+// //     nodeId2 = id_tri[1];
+// //     nodeId3 = id_tri[2];
+// 
+//     if ( nbTriangle > 1 ) {
+//       // get the nodes indices
+//       const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
+//       const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
+//       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
+//       try {
+//         OCC_CATCH_SIGNALS;
+//         tabID[i] = findShapeID( theMesh, n1, n2, n3, toMeshHoles );
+//         // -- 0020330: Pb with ghs3d as a submesh
+//         // check that found shape is to be meshed
+//         if ( tabID[i] > 0 ) {
+//           const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] );
+//           bool isToBeMeshed = false;
+//           for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS )
+//             isToBeMeshed = foundShape.IsSame( tabShape[ iS ]);
+//           if ( !isToBeMeshed )
+//             tabID[i] = HOLE_ID;
+//         }
+//         // END -- 0020330: Pb with ghs3d as a submesh
+// #ifdef _DEBUG_
+//         std::cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << std::endl;
+// #endif
+//       }
+//       catch ( Standard_Failure & ex)
+//       {
+// #ifdef _DEBUG_
+//         std::cout << i+1 << " subdomain: Exception caugt: " << ex.GetMessageString() << std::endl;
+// #endif
+//       }
+//       catch (...) {
+// #ifdef _DEBUG_
+//         std::cout << i+1 << " subdomain: unknown exception caught " << std::endl;
+// #endif
+//       }
+//     }
+//   }
+//   
+//   // ===========================
+//   // Fill the tabID array: END
+//   // ===========================
+//   
+// 
+//   tabRef[GmfVertices]       = 3;
+//   tabRef[GmfCorners]        = 1;
+//   tabRef[GmfEdges]          = 2;
+//   tabRef[GmfRidges]         = 1;
+//   tabRef[GmfTriangles]      = 3;
+// //   tabRef[GmfQuadrilaterals] = 4;
+//   tabRef[GmfTetrahedra]     = 4;
+// //   tabRef[GmfHexahedra]      = 8;
+//   
+//   SMDS_NodeIteratorPtr itOnGMFInputNode = theMeshDS->nodesIterator();
+//   while ( itOnGMFInputNode->more() )
+//     theMeshDS->RemoveNode( itOnGMFInputNode->next() );
+// 
+//   
+//   int nbVertices = GmfStatKwd(InpMsh, GmfVertices);
+//   int nbCorners = max(countShape( theMeshDS, TopAbs_VERTEX ) , GmfStatKwd(InpMsh, GmfCorners));
+//   int nbShapeEdge = countShape( theMeshDS, TopAbs_EDGE );
+// 
+//   tabCorner       = new TopoDS_Shape[ nbCorners ];
+//   tabEdge         = new TopoDS_Shape[ nbShapeEdge ];
+//   nodeAssigne     = new int[ nbVertices + 1 ];
+//   GMFNodeAssigne  = new int[ nbVertices + 1 ];
+//   GMFNode         = new SMDS_MeshNode*[ nbVertices + 1 ];
+// 
+//   getShape(theMeshDS, TopAbs_VERTEX, tabCorner);
+//   getShape(theMeshDS, TopAbs_EDGE,   tabEdge);
+// 
+//   std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
+//   for ( ; it != tabRef.end() ; ++it)
+//   {
+// //     int dummy;
+//     GmfKwdCod token = it->first;
+//     nbRef    = it->second;
+// 
+//     nbElem = GmfStatKwd(InpMsh, token);
+//     if (nbElem > 0) {
+//       GmfGotoKwd(InpMsh, token);
+//       std::cout << "Read " << nbElem;
+//     }
+//     else
+//       continue;
+// 
+//     int id[nbElem*tabRef[token]];
+//     int ghs3dShapeID[nbElem];
+// 
+//     if (token == GmfVertices) {
+//       std::cout << " vertices" << std::endl;
+//       int aGMFID;
+// 
+//       float VerTab_f[nbElem][3];
+//       double VerTab_d[nbElem][3];
+//       SMDS_MeshNode * aGMFNode;
+// 
+//       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
+//         aGMFID = iElem + 1;
+//         if (ver == GmfFloat) {
+//           GmfGetLin(InpMsh, token, &VerTab_f[nbElem][0], &VerTab_f[nbElem][1], &VerTab_f[nbElem][2], &ghs3dShapeID[iElem]);
+//           aGMFNode = theMeshDS->AddNode(VerTab_f[nbElem][0], VerTab_f[nbElem][1], VerTab_f[nbElem][2]);
+//         }
+//         else {
+//           GmfGetLin(InpMsh, token, &VerTab_d[nbElem][0], &VerTab_d[nbElem][1], &VerTab_d[nbElem][2], &ghs3dShapeID[iElem]);
+//           aGMFNode = theMeshDS->AddNode(VerTab_d[nbElem][0], VerTab_d[nbElem][1], VerTab_d[nbElem][2]);
+//         }
+//         GMFNode[ aGMFID ] = aGMFNode;
+//         nodeAssigne[ aGMFID ] = 0;
+//         GMFNodeAssigne[ aGMFID ] = 0;
+//       }
+//     }
+//     else if (token == GmfCorners && nbElem > 0) {
+//       std::cout << " corners" << std::endl;
+//       for ( int iElem = 0; iElem < nbElem; iElem++ )
+//         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
+//     }
+//     else if (token == GmfRidges && nbElem > 0) {
+//       std::cout << " ridges" << std::endl;
+//       for ( int iElem = 0; iElem < nbElem; iElem++ )
+//         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
+//     }
+//     else if (token == GmfEdges && nbElem > 0) {
+//       std::cout << " edges" << std::endl;
+//       for ( int iElem = 0; iElem < nbElem; iElem++ )
+//         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &ghs3dShapeID[iElem]);
+//     }
+//     else if (token == GmfTriangles && nbElem > 0) {
+//       std::cout << " triangles" << std::endl;
+//       for ( int iElem = 0; iElem < nbElem; iElem++ )
+//         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &ghs3dShapeID[iElem]);
+//     }
+// //     else if (token == GmfQuadrilaterals && nbElem > 0) {
+// //       std::cout << " Quadrilaterals" << std::endl;
+// //       for ( int iElem = 0; iElem < nbElem; iElem++ )
+// //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &ghs3dShapeID[iElem]);
+// //     }
+//     else if (token == GmfTetrahedra && nbElem > 0) {
+//       std::cout << " Tetrahedra" << std::endl;
+//       for ( int iElem = 0; iElem < nbElem; iElem++ )
+//         GmfGetLin(InpMsh, token, 
+//                   &id[iElem*tabRef[token]], 
+//                   &id[iElem*tabRef[token]+1], 
+//                   &id[iElem*tabRef[token]+2], 
+//                   &id[iElem*tabRef[token]+3], 
+//                   &ghs3dShapeID[iElem]);
+//     }
+// //     else if (token == GmfHexahedra && nbElem > 0) {
+// //       std::cout << " Hexahedra" << std::endl;
+// //       for ( int iElem = 0; iElem < nbElem; iElem++ )
+// //         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
+// //                   &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &ghs3dShapeID[iElem]);
+// //     }
+// 
+//     switch (token) {
+//     case GmfCorners:
+//     case GmfRidges:
+//     case GmfEdges:
+//     case GmfTriangles:
+// //     case GmfQuadrilaterals:
+//     case GmfTetrahedra:
+// //     case GmfHexahedra:
+//     {
+//       int nodeDim, shapeID, *nodeID;
+//       SMDS_MeshNode** node;
+// //       std::vector< SMDS_MeshNode* > enfNode( nbRef );
+//       SMDS_MeshElement * aGMFElement;
+//       
+//       node    = new SMDS_MeshNode*[nbRef];
+//       nodeID  = new int[ nbRef ];
+// 
+//       for ( int iElem = 0; iElem < nbElem; iElem++ )
+//       {
+//         for ( int iRef = 0; iRef < nbRef; iRef++ )
+//         {
+//           aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
+//           node  [ iRef ] = GMFNode[ aGMFNodeID ];
+//           nodeID[ iRef ] = aGMFNodeID;
+//         }
+// 
+//         switch (token)
+//         {
+//         case GmfCorners: {
+//           nodeDim = 1;
+//           gp_Pnt GMFPnt ( node[0]->X(), node[0]->Y(), node[0]->Z() );
+//           for ( int i=0; i<nbElem; i++ ) {
+//             aVertex = TopoDS::Vertex( tabCorner[i] );
+//             gp_Pnt aPnt = BRep_Tool::Pnt( aVertex );
+//             if ( aPnt.Distance( GMFPnt ) < epsilon )
+//               break;
+//           }
+//           break;
+//         }
+//         case GmfEdges: {
+//           nodeDim = 2;
+//           aGMFElement = theMeshDS->AddEdge( node[0], node[1] );
+//           int iNode = 1;
+//           if ( GMFNodeAssigne[ nodeID[0] ] == 0 || GMFNodeAssigne[ nodeID[0] ] == 2 )
+//             iNode = 0;
+//           shapeID = findEdgeID( node[iNode], theMeshDS, nbShapeEdge, tabEdge );
+//           break;
+//         }
+//         case GmfRidges:
+//           break;
+//         case GmfTriangles: {
+//           nodeDim = 3;
+//           aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2]);
+//           shapeID = -1;
+//           break;
+//         }
+// //         case GmfQuadrilaterals: {
+// //           nodeDim = 4;
+// //           aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2], node[3] );
+// //           shapeID = -1;
+// //           break;
+// //         }
+//         case GmfTetrahedra: {
+//           
+//           // IN WORK
+//           TopoDS_Shape aSolid;
+//           // We always run GHS3D with "to mesh holes"==TRUE but we must not create
+//           // tetras within holes depending on hypo option,
+//           // so we first check if aTet is inside a hole and then create it 
+//           if ( nbTriangle > 1 ) {
+//             tetraShapeID = HOLE_ID; // negative tetraShapeID means not to create tetras if !toMeshHoles
+//             int aGhs3dShapeID = ghs3dShapeID[iElem] - IdShapeRef;
+//             if ( tabID[ aGhs3dShapeID ] == 0 ) {
+//               TopAbs_State state;
+//               aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
+//               if ( toMeshHoles || state == TopAbs_IN )
+//                 tetraShapeID = theMeshDS->ShapeToIndex( aSolid );
+//               tabID[ aGhs3dShapeID ] = tetraShapeID;
+//             }
+//             else
+//               tetraShapeID = tabID[ aGhs3dShapeID ];
+//           }
+//           else if ( nbShape > 1 ) {
+//             // Case where nbTriangle == 1 while nbShape == 2 encountered
+//             // with compound of 2 boxes and "To mesh holes"==False,
+//             // so there are no subdomains specified for each tetrahedron.
+//             // Try to guess a solid by a node already bound to shape
+//             tetraShapeID = 0;
+//             for ( int i=0; i<4 && tetraShapeID==0; i++ ) {
+//               if ( nodeAssigne[ nodeID[i] ] == 1 &&
+//                   node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
+//                   node[i]->getshapeId() > 1 )
+//               {
+//                 tetraShapeID = node[i]->getshapeId();
+//               }
+//             }
+//             if ( tetraShapeID==0 ) {
+//               aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
+//               tetraShapeID = theMeshDS->ShapeToIndex( aSolid );
+//             }
+//           }
+//           // set new nodes and tetrahedron onto the shape
+//           for ( int i=0; i<4; i++ ) {
+//             if ( nodeAssigne[ nodeID[i] ] == 0 ) {
+//               if ( tetraShapeID != HOLE_ID )
+//                 theMeshDS->SetNodeInVolume( node[i], tetraShapeID );
+//               nodeAssigne[ nodeID[i] ] = tetraShapeID;
+//             }
+//           }
+//           if ( toMeshHoles || tetraShapeID != HOLE_ID ) {
+//             aGMFElement = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
+//             theMeshDS->SetMeshElementOnShape( aGMFElement, tetraShapeID );
+//           }
+//           
+//           // IN WORK
+//           
+//           nodeDim = 5;
+//           break;
+//         }
+// //         case GmfHexahedra: {
+// //           nodeDim = 6;
+// //           aGMFElement = theMeshDS->AddVolume( node[0], node[3], node[2], node[1],
+// //                                             node[4], node[7], node[6], node[5] );
+// //           break;
+// //         }
+//         default: continue;
+//         }
+//         if (token != GmfRidges)
+//         {
+//           for ( int i=0; i<nbRef; i++ ) {
+//               if ( GMFNodeAssigne[ nodeID[i] ] == 0 ) {
+//                 if      ( token == GmfCorners )   theMeshDS->SetNodeOnVertex( node[0], aVertex );
+//                 else if ( token == GmfEdges )     theMeshDS->SetNodeOnEdge( node[i], shapeID );
+//                 else if ( token == GmfTriangles ) theMeshDS->SetNodeOnFace( node[i], shapeID );
+//                 GMFNodeAssigne[ nodeID[i] ] = nodeDim;
+//               }
+//             }
+//             if ( token != "Corners" )
+//               theMeshDS->SetMeshElementOnShape( aGMFElement, shapeID );
+//         }
+//       } // for
+//       
+//       if ( !toMeshHoles ) {
+//         map <int,const SMDS_MeshNode*>::iterator itOnNode = theGhs3dIdToNodeMap.find( nbVertices-(nbEnforcedVertices+nbEnforcedNodes) );
+//         for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
+//           if ( nodeAssigne[ itOnNode->first ] == HOLE_ID )
+//             theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
+//         }
+//       }
+//       
+//       delete [] node;
+//       delete [] nodeID;
+//       break;
+//       } // case GmfTetrahedra
+//     } // switch(token)
+//   } // for
+//   cout << std::endl;
+//   
+// #ifdef WNT
+//   UnmapViewOfFile(mapPtr);
+//   CloseHandle(hMapObject);
+//   CloseHandle(fd);
+// #else
+//   munmap(mapPtr, length);
+// #endif
+//   close(fileOpen);
+//   
+//   delete [] tabID;
+//   delete [] tabCorner;
+//   delete [] tabEdge;
+//   delete [] nodeAssigne;
+//   delete [] GMFNodeAssigne;
+//   delete [] GMFNode;
+//   
+//   return true;
+// }
+
+
+//=======================================================================
+//function : readGMFFile
+//purpose  : read GMF file w/o geometry associated to mesh
+//=======================================================================
+
+
+static bool readGMFFile(const char* theFile,
+#ifdef WITH_SMESH_CANCEL_COMPUTE
+                        GHS3DPlugin_GHS3D*  theAlgo,
+#endif 
+                        SMESH_MesherHelper* theHelper,
+                        TIDSortedNodeSet &  theEnforcedNodes,
+                        TIDSortedElemSet &  theEnforcedTriangles,
+                        TIDSortedElemSet &  theEnforcedQuadrangles)
+{
+  SMESHDS_Mesh* theMesh = theHelper->GetMeshDS();
+
+  // ---------------------------------
+  // Read generated elements and nodes
+  // ---------------------------------
+
+  int nbElem = 0, nbRef = 0;
+  int aGMFNodeID = 0, shapeID;
+  int *nodeAssigne;
+  SMDS_MeshNode** GMFNode;
+  std::map <GmfKwdCod,int> tabRef;
+
+  tabRef[GmfVertices]       = 3;
+  tabRef[GmfCorners]        = 1;
+  tabRef[GmfEdges]          = 2;
+  tabRef[GmfRidges]         = 1;
+  tabRef[GmfTriangles]      = 3;
+  tabRef[GmfQuadrilaterals] = 4;
+  tabRef[GmfTetrahedra]     = 4;
+  tabRef[GmfHexahedra]      = 8;
+
+  theHelper->GetMesh()->Clear();
+
+  int ver, dim;
+  MESSAGE("Read " << theFile << " file");
+  int InpMsh = GmfOpenMesh(theFile, GmfRead, &ver, &dim);
+  if (!InpMsh)
+    return false;
+
+  int nbVertices = GmfStatKwd(InpMsh, GmfVertices);
+  GMFNode = new SMDS_MeshNode*[ nbVertices + 1 ];
+  nodeAssigne = new int[ nbVertices + 1 ];
+
+  std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
+  for ( ; it != tabRef.end() ; ++it)
+  {
+#ifdef WITH_SMESH_CANCEL_COMPUTE
+    if(theAlgo->computeCanceled()) {
+      GmfCloseMesh(InpMsh);
+      delete [] GMFNode;
+      delete [] nodeAssigne;
+      return false;
+    }
+#endif
+    int dummy;
+    GmfKwdCod token = it->first;
+    nbRef    = it->second;
+
+    nbElem = GmfStatKwd(InpMsh, token);
+    if (nbElem > 0) {
+      GmfGotoKwd(InpMsh, token);
+      std::cout << "Read " << nbElem;
+    }
+    else
+      continue;
+
+    int id[nbElem*tabRef[token]];
+
+    if (token == GmfVertices) {
+      std::cout << " vertices" << std::endl;
+      int aGMFID;
+
+      float VerTab_f[nbElem][3];
+      double VerTab_d[nbElem][3];
+      SMDS_MeshNode * aGMFNode;
+
+      for ( int iElem = 0; iElem < nbElem; iElem++ ) {
+#ifdef WITH_SMESH_CANCEL_COMPUTE
+        if(theAlgo->computeCanceled()) {
+          GmfCloseMesh(InpMsh);
+          delete [] GMFNode;
+          delete [] nodeAssigne;
+          return false;
+        }
+#endif
+        aGMFID = iElem + 1;
+        if (ver == GmfFloat) {
+          GmfGetLin(InpMsh, token, &VerTab_f[nbElem][0], &VerTab_f[nbElem][1], &VerTab_f[nbElem][2], &dummy);
+          aGMFNode = theMesh->AddNode(VerTab_f[nbElem][0], VerTab_f[nbElem][1], VerTab_f[nbElem][2]);
+        }
+        else {
+          GmfGetLin(InpMsh, token, &VerTab_d[nbElem][0], &VerTab_d[nbElem][1], &VerTab_d[nbElem][2], &dummy);
+          aGMFNode = theMesh->AddNode(VerTab_d[nbElem][0], VerTab_d[nbElem][1], VerTab_d[nbElem][2]);
+        }
+        GMFNode[ aGMFID ] = aGMFNode;
+        nodeAssigne[ aGMFID ] = 0;
+      }
+    }
+    else if (token == GmfCorners && nbElem > 0) {
+      std::cout << " corners" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
+    }
+    else if (token == GmfRidges && nbElem > 0) {
+      std::cout << " ridges" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
+    }
+    else if (token == GmfEdges && nbElem > 0) {
+      std::cout << " edges" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &dummy);
+    }
+    else if (token == GmfTriangles && nbElem > 0) {
+      std::cout << " triangles" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &dummy);
+    }
+    else if (token == GmfQuadrilaterals && nbElem > 0) {
+      std::cout << " Quadrilaterals" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &dummy);
+    }
+    else if (token == GmfTetrahedra && nbElem > 0) {
+      std::cout << " Tetrahedra" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &dummy);
+    }
+    else if (token == GmfHexahedra && nbElem > 0) {
+      std::cout << " Hexahedra" << std::endl;
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+        GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
+                  &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &dummy);
+    }
+    std::cout << std::endl;
+
+    switch (token) {
+    case GmfCorners:
+    case GmfRidges:
+    case GmfEdges:
+    case GmfTriangles:
+    case GmfQuadrilaterals:
+    case GmfTetrahedra:
+    case GmfHexahedra: {
+      std::vector< SMDS_MeshNode* > node( nbRef );
+      std::vector< int >          nodeID( nbRef );
+      std::vector< SMDS_MeshNode* > enfNode( nbRef );
+
+      for ( int iElem = 0; iElem < nbElem; iElem++ )
+      {
+#ifdef WITH_SMESH_CANCEL_COMPUTE
+        if(theAlgo->computeCanceled()) {
+          GmfCloseMesh(InpMsh);
+          delete [] GMFNode;
+          delete [] nodeAssigne;
+          return false;
+        }
+#endif
+        for ( int iRef = 0; iRef < nbRef; iRef++ )
+        {
+          aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
+          node  [ iRef ] = GMFNode[ aGMFNodeID ];
+          nodeID[ iRef ] = aGMFNodeID;
+        }
+
+        switch (token)
+        {
+        case GmfEdges:
+          theHelper->AddEdge( node[0], node[1] ); break;
+        case GmfTriangles: {
+          theMesh->AddFace( node[0], node[1], node[2]);
+          break;
+        }
+        case GmfQuadrilaterals: {
+          theMesh->AddFace( node[0], node[1], node[2], node[3] );
+          break;
+        }
+        case GmfTetrahedra:
+          theHelper->AddVolume( node[0], node[1], node[2], node[3] ); break;
+        case GmfHexahedra:
+          theHelper->AddVolume( node[0], node[3], node[2], node[1],
+                                node[4], node[7], node[6], node[5] ); break;
+        default: continue;
+        }
+        if ( token == GmfTriangles || token == GmfQuadrilaterals ) // "Quadrilaterals" and "Triangles"
+          for ( int iRef = 0; iRef < nbRef; iRef++ )
+            nodeAssigne[ nodeID[ iRef ]] = 1;
+      }
+      break;
+    }
+    }
+  }
+
+  shapeID = theHelper->GetSubShapeID();
+  for ( int i = 0; i < nbVertices; ++i ) {
+#ifdef WITH_SMESH_CANCEL_COMPUTE
+    if(theAlgo->computeCanceled()) {
+      GmfCloseMesh(InpMsh);
+      delete [] GMFNode;
+      delete [] nodeAssigne;
+      return false;
+    }
+#endif
+    if ( !nodeAssigne[ i+1 ])
+      theMesh->SetNodeInVolume( GMFNode[ i+1 ], shapeID );
+  }
+
+  GmfCloseMesh(InpMsh);
+  delete [] GMFNode;
+  delete [] nodeAssigne;
+  return true;
+}
+
+static bool writeGMFFile(const char*   theMeshFileName,
+                         const char*   theRequiredFileName,
+                         const char*   theSolFileName,
+                         const SMESH_ProxyMesh&           theProxyMesh,
+                         SMESH_Mesh *                     theMesh,
+                         vector <const SMDS_MeshNode*> &  theNodeByGhs3dId,
+                         vector <const SMDS_MeshNode*> &  theEnforcedNodeByGhs3dId,
+                         TIDSortedNodeSet &               theEnforcedNodes,
+                         TIDSortedElemSet &               theEnforcedEdges,
+                         TIDSortedElemSet &               theEnforcedTriangles,
+                         TIDSortedElemSet &               theEnforcedQuadrangles,
+                         GHS3DPlugin_Hypothesis::TEnforcedVertexValues & theEnforcedVertices)
+{
+  MESSAGE("writeGMFFile w/o geometry");
+  int idx, idxRequired, idxSol;
+  const int dummyint = 0;
+  GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
+  std::vector<double> enfVertexSizes;
+  const SMDS_MeshElement* elem;
+  TIDSortedElemSet anElemSet, anEnforcedEdgeSet, anEnforcedTriangleSet, anEnforcedQuadrangleSet;
+  SMDS_ElemIteratorPtr nodeIt;
+  map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap, anEnforcedNodeToGhs3dIdMap;
+  TIDSortedElemSet::iterator elemIt;
+  bool isOK;
+  auto_ptr< SMESH_ElementSearcher > pntCls ( SMESH_MeshEditor( theMesh ).GetElementSearcher());
+  
+  int nbEnforcedVertices = theEnforcedVertices.size();
+//  int nbEnforcedNodes    = theEnforcedNodes.size();
+//  int nbEnforcedEdges       = theEnforcedEdges.size();
+//  int nbEnforcedTriangles   = theEnforcedTriangles.size();
+//  int nbEnforcedQuadrangles = theEnforcedQuadrangles.size();
+  
+  // count faces
+  int nbFaces = theProxyMesh.NbFaces();
+  int nbNodes;
+
+  if ( nbFaces == 0 )
+    return false;
+  
+  idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+  if (!idx)
+    return false;
+  
+  /* ========================== FACES ========================== */
+  /* TRIANGLES ========================== */
+  SMDS_ElemIteratorPtr eIt = theProxyMesh.GetFaces();
+  while ( eIt->more() )
+  {
+    elem = eIt->next();
+    anElemSet.insert(elem);
+    // NODE_NB_1 NODE_NB_2 ...
+    nodeIt = elem->nodesIterator();
+    nbNodes = elem->NbCornerNodes();
+    while ( nodeIt->more() && nbNodes-- )
+    {
+      // find GHS3D ID
+      const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+      int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
+      aNodeToGhs3dIdMap.insert( make_pair( node, newId ));
+    }
+  }
+  
+  /* EDGES ========================== */
+  
+  // Iterate over the enforced edges
+  for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
+    elem = (*elemIt);
+    isOK = true;
+    nodeIt = elem->nodesIterator();
+    nbNodes = elem->NbCornerNodes();
+    while ( nodeIt->more() && nbNodes-- ) {
+      // find GHS3D ID
+      const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+      // Test if point is inside shape to mesh
+      gp_Pnt myPoint(node->X(),node->Y(),node->Z());
+      TopAbs_State result = pntCls->GetPointState( myPoint );
+      if ( result != TopAbs_IN ) {
+        isOK = false;
+        break;
+      }
+    }
+    if (isOK) {
+      nodeIt = elem->nodesIterator();
+      nbNodes = elem->NbCornerNodes();
+      while ( nodeIt->more() && nbNodes-- ) {
+        // find GHS3D ID
+        const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+        int newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
+        anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
+      }
+      anEnforcedEdgeSet.insert(elem);
+    }
+  }
+  
+  /* ENFORCED TRIANGLES ========================== */
+  
+  // Iterate over the enforced triangles
+  for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
+    elem = (*elemIt);
+    isOK = true;
+    nodeIt = elem->nodesIterator();
+    nbNodes = elem->NbCornerNodes();
+    while ( nodeIt->more() && nbNodes-- ) {
+      // find GHS3D ID
+      const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+      // Test if point is inside shape to mesh
+      gp_Pnt myPoint(node->X(),node->Y(),node->Z());
+      TopAbs_State result = pntCls->GetPointState( myPoint );
+      if ( result != TopAbs_IN ) {
+        isOK = false;
+        break;
+      }
+    }
+    if (isOK) {
+      nodeIt = elem->nodesIterator();
+      nbNodes = elem->NbCornerNodes();
+      while ( nodeIt->more() && nbNodes-- ) {
+        // find GHS3D ID
+        const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+        int newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
+        anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
+      }
+      anEnforcedTriangleSet.insert(elem);
+    }
+  }
+  
+  /* ENFORCED QUADRANGLES ========================== */
+  
+    // Iterate over the enforced quadrangles
+  for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) {
+    elem = (*elemIt);
+    isOK = true;
+    nodeIt = elem->nodesIterator();
+    nbNodes = elem->NbCornerNodes();
+    while ( nodeIt->more() && nbNodes-- ) {
+      // find GHS3D ID
+      const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+      // Test if point is inside shape to mesh
+      gp_Pnt myPoint(node->X(),node->Y(),node->Z());
+      TopAbs_State result = pntCls->GetPointState( myPoint );
+      if ( result != TopAbs_IN ) {
+        isOK = false;
+        break;
+      }
+    }
+    if (isOK) {
+      nodeIt = elem->nodesIterator();
+      nbNodes = elem->NbCornerNodes();
+      while ( nodeIt->more() && nbNodes-- ) {
+        // find GHS3D ID
+        const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+        int newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
+        anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
+      }
+      anEnforcedQuadrangleSet.insert(elem);
+    }
+  }
+  
+  
+  // put nodes to theNodeByGhs3dId vector
+  std::cout << "aNodeToGhs3dIdMap.size(): "<<aNodeToGhs3dIdMap.size()<<std::endl;
+  theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
+  map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
+  for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
+  {
+//     std::cout << "n2id->first: "<<n2id->first<<std::endl;
+    theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
+  }
+
+  // put nodes to anEnforcedNodeToGhs3dIdMap vector
+  std::cout << "anEnforcedNodeToGhs3dIdMap.size(): "<<anEnforcedNodeToGhs3dIdMap.size()<<std::endl;
+  theEnforcedNodeByGhs3dId.resize( anEnforcedNodeToGhs3dIdMap.size() );
+  n2id = anEnforcedNodeToGhs3dIdMap.begin();
+  for ( ; n2id != anEnforcedNodeToGhs3dIdMap.end(); ++ n2id)
+  {
+//     std::cout << "n2id->first: "<<n2id->first<<std::endl;
+    theEnforcedNodeByGhs3dId[ n2id->second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // ghs3d ids count from 1
+  }
+  
+  
+  /* ========================== NODES ========================== */
+  vector<const SMDS_MeshNode*> theOrderedNodes, theRequiredNodes;
+  std::set< std::vector<double> > nodesCoords;
+  vector<const SMDS_MeshNode*>::const_iterator ghs3dNodeIt = theNodeByGhs3dId.begin();
+  vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
+  
+  std::cout << theNodeByGhs3dId.size() << " nodes from mesh ..." << std::endl;
+  for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
+  {
+    const SMDS_MeshNode* node = *ghs3dNodeIt;
+    std::vector<double> coords;
+    coords.push_back(node->X());
+    coords.push_back(node->Y());
+    coords.push_back(node->Z());
+    nodesCoords.insert(coords);
+    theOrderedNodes.push_back(node);
+  }
+  
+  // Iterate over the enforced nodes
+  TIDSortedNodeSet::const_iterator enfNodeIt;
+  std::cout << theEnforcedNodes.size() << " nodes from enforced nodes ..." << std::endl;
+  for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt)
+  {
+    const SMDS_MeshNode* node = *enfNodeIt;
+    std::vector<double> coords;
+    coords.push_back(node->X());
+    coords.push_back(node->Y());
+    coords.push_back(node->Z());
+    
+    if (nodesCoords.find(coords) != nodesCoords.end()) {
+      std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z() << " found" << std::endl;
+      continue;
+    }
+
+    if (theEnforcedVertices.find(coords) != theEnforcedVertices.end())
+      continue;
+    
+    // Test if point is inside shape to mesh
+    gp_Pnt myPoint(node->X(),node->Y(),node->Z());
+    TopAbs_State result = pntCls->GetPointState( myPoint );
+    if ( result != TopAbs_IN )
+      continue;
+    
+    nodesCoords.insert(coords);
+    theOrderedNodes.push_back(node);
+    theRequiredNodes.push_back(node);
+  }
+  // Iterate over the enforced nodes given by enforced elements
+  ghs3dNodeIt = theEnforcedNodeByGhs3dId.begin();
+  after  = theEnforcedNodeByGhs3dId.end();
+  std::cout << theEnforcedNodeByGhs3dId.size() << " nodes from enforced elements ..." << std::endl;
+  for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
+  {
+    const SMDS_MeshNode* node = *ghs3dNodeIt;
+    std::vector<double> coords;
+    coords.push_back(node->X());
+    coords.push_back(node->Y());
+    coords.push_back(node->Z());
+    
+    if (nodesCoords.find(coords) != nodesCoords.end()) {
+      std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z() << " found" << std::endl;
+      continue;
+    }
+    
+    if (theEnforcedVertices.find(coords) != theEnforcedVertices.end())
+      continue;
+    
+    nodesCoords.insert(coords);
+    theOrderedNodes.push_back(node);
+    theRequiredNodes.push_back(node);
+  }
+  
+  int requiredNodes = theRequiredNodes.size();
+  int solSize = 0;
+  std::vector<std::vector<double> > ReqVerTab;
+  if (nbEnforcedVertices) {
+//    ReqVerTab.clear();
+    std::cout << theEnforcedVertices.size() << " nodes from enforced vertices ..." << std::endl;
+    // Iterate over the enforced vertices
+    for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
+      double x = vertexIt->first[0];
+      double y = vertexIt->first[1];
+      double z = vertexIt->first[2];
+      // Test if point is inside shape to mesh
+      gp_Pnt myPoint(x,y,z);
+      TopAbs_State result = pntCls->GetPointState( myPoint );
+      if ( result != TopAbs_IN )
+        continue;
+      std::vector<double> coords;
+      coords.push_back(x);
+      coords.push_back(y);
+      coords.push_back(z);
+      ReqVerTab.push_back(coords);
+      enfVertexSizes.push_back(vertexIt->second);
+      solSize++;
+    }
+  }
+
+  // GmfVertices
+  std::cout << "Begin writting required nodes in GmfVertices" << std::endl;
+  GmfSetKwd(idx, GmfVertices, theOrderedNodes.size()+solSize);
+  for (ghs3dNodeIt = theOrderedNodes.begin();ghs3dNodeIt != theOrderedNodes.end();++ghs3dNodeIt)
+    GmfSetLin(idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint);
+  for (int i=0;i<solSize;i++) {
+    std::cout << ReqVerTab[i][0] <<" "<< ReqVerTab[i][1] << " "<< ReqVerTab[i][2] << std::endl;
+    GmfSetLin(idx, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
+  }
+  std::cout << "End writting required nodes in GmfVertices" << std::endl;
+
+  if (requiredNodes + solSize) {
+    GmfSetKwd(idx, GmfRequiredVertices, requiredNodes+solSize);
+    if (requiredNodes) {
+      std::cout << "Begin writting required nodes in GmfRequiredVertices" << std::endl;
+      int startId = theOrderedNodes.size()-requiredNodes+1;
+      std::cout << "startId: " << startId << std::endl;
+      for (int i=0;i<requiredNodes;i++)
+        GmfSetLin(idx, GmfRequiredVertices, startId+i);
+      std::cout << "End writting required nodes in GmfRequiredVertices" << std::endl;
+    }
+    if (solSize) {
+      std::cout << "Begin writting required vertices in GmfRequiredVertices" << std::endl;
+      int startId = theOrderedNodes.size()+1;
+      std::cout << "startId: " << startId << std::endl;
+      for (int i=0;i<solSize;i++)
+        GmfSetLin(idx, GmfRequiredVertices, startId+i);
+      std::cout << "End writting required vertices in GmfRequiredVertices" << std::endl;
+
+      std::cout << "Begin writting in sol file" << std::endl;
+      idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+      if (!idxSol) {
+        GmfCloseMesh(idx);
+        if (idxRequired)
+          GmfCloseMesh(idxRequired);
+        return false;
+      }
+      int TypTab[] = {GmfSca};
+      GmfSetKwd(idxSol, GmfSolAtVertices, solSize, 1, TypTab);
+      for (int i=0;i<solSize;i++) {
+        std::cout << "enfVertexSizes.at(i): " << enfVertexSizes.at(i) << std::endl;
+        double solTab[] = {enfVertexSizes.at(i)};
+        GmfSetLin(idxSol, GmfSolAtVertices, solTab);
+      }
+      std::cout << "End writting in sol file" << std::endl;
+    }
+  }
+
+//  // GmfRequiredVertices + GmfSolAtVertices
+////  std::cout << "theRequiredNodes.size() + solSize: " << theRequiredNodes.size()+ solSize << std::endl;
+////  std::cout << "theRequiredNodes.size(): " << theRequiredNodes.size() << std::endl;
+//  std::cout << "solSize: " << solSize << std::endl;
+////  if (theRequiredNodes.size()+ solSize) {
+////    GmfSetKwd(idx, GmfRequiredVertices, theRequiredNodes.size()+solSize);
+////
+////    if (theRequiredNodes.size()) {
+////      std::cout << "Begin writting required nodes in GmfRequiredVertices" << std::endl;
+////      int startId = theOrderedNodes.size()-theRequiredNodes.size();
+////      std::cout << "startId: " << startId << std::endl;
+////      for (int i=1;i<=theRequiredNodes.size();i++)
+////        GmfSetLin(idx, GmfRequiredVertices, startId+i);
+////      std::cout << "End writting required nodes in GmfRequiredVertices" << std::endl;
+////    }
+//
+//    if (solSize) {
+//      std::cout << "Begin writting in sol file" << std::endl;
+//      GmfSetKwd(idx, GmfRequiredVertices, solSize);
+//      idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+//      if (!idxSol) {
+//        GmfCloseMesh(idx);
+//        if (idxRequired)
+//          GmfCloseMesh(idxRequired);
+//        return false;
+//      }
+//
+//      int TypTab[] = {GmfSca};
+////      GmfSetKwd(idxRequired, GmfVertices, solSize);
+//      GmfSetKwd(idxSol, GmfSolAtVertices, solSize, 1, TypTab);
+//
+//      for (int i=0;i<solSize;i++) {
+//        double solTab[] = {enfVertexSizes.at(i)};
+//        GmfSetLin(idx, GmfRequiredVertices, theOrderedNodes.size()+i+1);
+//        GmfSetLin(idxSol, GmfSolAtVertices, solTab);
+////      GmfSetLin(idxRequired, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
+//      }
+//      std::cout << "End writting in sol file" << std::endl;
+//    }
+//
+////  }
+
+  int nedge[2], ntri[3], nquad[4];
+  
+  int usedEnforcedEdges = 0;
+  if (anEnforcedEdgeSet.size()) {
+//    idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+//    if (!idxRequired)
+//      return false;
+    GmfSetKwd(idx, GmfEdges, anEnforcedEdgeSet.size());
+//    GmfSetKwd(idxRequired, GmfEdges, anEnforcedEdgeSet.size());
+    for(elemIt = anEnforcedEdgeSet.begin() ; elemIt != anEnforcedEdgeSet.end() ; ++elemIt) {
+      elem = (*elemIt);
+      nodeIt = elem->nodesIterator();
+      int index=0;
+      while ( nodeIt->more() ) {
+        // find GHS3D ID
+        const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+        map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
+        if (it == anEnforcedNodeToGhs3dIdMap.end())
+          throw "Node not found";
+        nedge[index] = it->second;
+        index++;
+      }
+      GmfSetLin(idx, GmfEdges, nedge[0], nedge[1], dummyint);
+//      GmfSetLin(idxRequired, GmfEdges, nedge[0], nedge[1], dummyint);
+      usedEnforcedEdges++;
+    }
+//    GmfCloseMesh(idxRequired);
+  }
+
+  if (usedEnforcedEdges) {
+    GmfSetKwd(idx, GmfRequiredEdges, usedEnforcedEdges);
+    for (int enfID=1;enfID<=usedEnforcedEdges;enfID++) {
+      GmfSetLin(idx, GmfRequiredEdges, enfID);
+    }
+  }
+
+  if (anElemSet.size()+anEnforcedTriangleSet.size()) {
+    GmfSetKwd(idx, GmfTriangles, anElemSet.size()+anEnforcedTriangleSet.size());
+    for(elemIt = anElemSet.begin() ; elemIt != anElemSet.end() ; ++elemIt) {
+      elem = (*elemIt);
+      nodeIt = elem->nodesIterator();
+      int index=0;
+      while ( nodeIt->more() ) {
+        // find GHS3D ID
+        const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+        map< const SMDS_MeshNode*,int >::iterator it = aNodeToGhs3dIdMap.find(node);
+        if (it == aNodeToGhs3dIdMap.end())
+          throw "Node not found";
+        ntri[index] = it->second;
+        index++;
+      }
+      GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint);
+    }
+    if (anEnforcedTriangleSet.size()) {
+      int usedEnforcedTriangles = 0;
+      for(elemIt = anEnforcedTriangleSet.begin() ; elemIt != anEnforcedTriangleSet.end() ; ++elemIt) {
+        elem = (*elemIt);
+        nodeIt = elem->nodesIterator();
+        int index=0;
+        while ( nodeIt->more() ) {
+          // find GHS3D ID
+          const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+          map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
+          if (it == anEnforcedNodeToGhs3dIdMap.end())
+            throw "Node not found";
+          ntri[index] = it->second;
+          index++;
+        }
+        GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint);
+        usedEnforcedTriangles++;
+      }
+      if (usedEnforcedTriangles) {
+        GmfSetKwd(idx, GmfRequiredTriangles, usedEnforcedTriangles);
+        for (int enfID=1;enfID<=usedEnforcedTriangles;enfID++)
+          GmfSetLin(idx, GmfRequiredTriangles, anElemSet.size()+enfID);
+      }
+    }
+  }
+
+  if (anEnforcedQuadrangleSet.size()) {
+    int usedEnforcedQuadrilaterals = 0;
+    GmfSetKwd(idx, GmfQuadrilaterals, anEnforcedQuadrangleSet.size());
+    for(elemIt = anEnforcedQuadrangleSet.begin() ; elemIt != anEnforcedQuadrangleSet.end() ; ++elemIt) {
+      elem = (*elemIt);
+      nodeIt = elem->nodesIterator();
+      int index=0;
+      while ( nodeIt->more() ) {
+        // find GHS3D ID
+        const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+        map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
+        if (it == anEnforcedNodeToGhs3dIdMap.end())
+          throw "Node not found";
+        nquad[index] = it->second;
+        index++;
+      }
+      GmfSetLin(idx, GmfQuadrilaterals, nquad[0], nquad[1], nquad[2], nquad[3], dummyint);
+      usedEnforcedQuadrilaterals++;
+    }
+    if (usedEnforcedQuadrilaterals) {
+      GmfSetKwd(idx, GmfRequiredQuadrilaterals, usedEnforcedQuadrilaterals);
+      for (int enfID=1;enfID<=usedEnforcedQuadrilaterals;enfID++)
+        GmfSetLin(idx, GmfRequiredQuadrilaterals, enfID);
+    }
+  }
+
+  GmfCloseMesh(idx);
+  if (idxRequired)
+    GmfCloseMesh(idxRequired);
+  if (idxSol)
+    GmfCloseMesh(idxSol);
+  
+  return true;
+  
+}
+
+static bool writeGMFFile(const char*   theMeshFileName,
+                         const char*   theRequiredFileName,
+                         const char*   theSolFileName,
+                         SMESH_MesherHelper&              theHelper,
+                         const SMESH_ProxyMesh&           theProxyMesh,
+                         map <int,int> &                  theSmdsToGhs3dIdMap,
+                         map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap,
+                         TIDSortedNodeSet &               theEnforcedNodes,
+                         TIDSortedElemSet &               theEnforcedEdges,
+                         TIDSortedElemSet &               theEnforcedTriangles,
+                         TIDSortedElemSet &               theEnforcedQuadrangles,
+                         GHS3DPlugin_Hypothesis::TEnforcedVertexValues & theEnforcedVertices)
+{
+  MESSAGE("writeGMFFile with geometry");
+  int idx, nbv, nbev, nben, aGhs3dID = 0;
+  const int dummyint = 0;
+  GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
+  std::vector<double> enfVertexSizes;
+  TIDSortedNodeSet::const_iterator enfNodeIt;
+  const SMDS_MeshNode* node;
+  SMDS_NodeIteratorPtr nodeIt;
+  std::map<int,int> theNodeId2NodeIndexMap;
+
+  idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+  if (!idx)
+    return false;
+  
+  SMESHDS_Mesh * theMeshDS = theHelper.GetMeshDS();
+  
+  /* ========================== NODES ========================== */
+  // NB_NODES
+  nbv = theMeshDS->NbNodes();
+  if ( nbv == 0 )
+    return false;
+  nbev = theEnforcedVertices.size();
+  nben = theEnforcedNodes.size();
+  
+  nodeIt = theMeshDS->nodesIterator();
+  
+  // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
+  // The problem is in nodes on degenerated edges, we need to skip them
+  if ( theHelper.HasDegeneratedEdges() )
+  {
+    // here we decrease total nb of nodes by nb of nodes on degenerated edges
+    set<int> checkedSM;
+    for (TopExp_Explorer e(theMeshDS->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
+    {
+      SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
+      if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() )) {
+        if ( sm->GetSubMeshDS() )
+          nbv -= sm->GetSubMeshDS()->NbNodes();
+      }
+    }
+  }
+  
+  const bool isQuadMesh = 
+    theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
+    theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
+    theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC );
+  if ( isQuadMesh )
+  {
+    // descrease nbNodes by nb of medium nodes
+    while ( nodeIt->more() )
+    {
+      node = nodeIt->next();
+      if ( !theHelper.IsDegenShape( node->getshapeId() ))
+        nbv -= int( theHelper.IsMedium( node ));
+    }
+    nodeIt = theMeshDS->nodesIterator();
+  }
+  
+  std::vector<std::vector<double> > VerTab;
+  VerTab.clear();
+  std::vector<double> aVerTab;
+  // Loop from 1 to NB_NODES
+
+  while ( nodeIt->more() )
+  {
+    node = nodeIt->next();
+    if (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238
+        theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674
+      continue;
+
+    aVerTab.clear();
+    aVerTab.push_back(node->X());
+    aVerTab.push_back(node->Y());
+    aVerTab.push_back(node->Z());
+    VerTab.push_back(aVerTab);
+    aGhs3dID++;
+    theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
+    theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
+  }
+  
+  /* ENFORCED NODES ========================== */
+  if (nben) {
+    for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt) {
+      double x = (*enfNodeIt)->X();
+      double y = (*enfNodeIt)->Y();
+      double z = (*enfNodeIt)->Z();
+      // Test if point is inside shape to mesh
+      gp_Pnt myPoint(x,y,z);
+      BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh());
+      scl.Perform(myPoint, 1e-7);
+      TopAbs_State result = scl.State();
+      if ( result != TopAbs_IN )
+        continue;
+      std::vector<double> coords;
+      coords.push_back(x);
+      coords.push_back(y);
+      coords.push_back(z);
+      if (theEnforcedVertices.find(coords) != theEnforcedVertices.end())
+        continue;
+      aVerTab.clear();
+      aVerTab.push_back(x);
+      aVerTab.push_back(y);
+      aVerTab.push_back(z);
+      VerTab.push_back(aVerTab);
+      aGhs3dID++;
+      theNodeId2NodeIndexMap.insert( make_pair( (*enfNodeIt)->GetID(), aGhs3dID ));
+    }
+  }
+  
+  /* Write vertices number */
+  MESSAGE("Number of vertices: "<<aGhs3dID);
+  MESSAGE("Size of vector: "<<VerTab.size());
+  GmfSetKwd(idx, GmfVertices, aGhs3dID);
+  for (int i=0;i<aGhs3dID;i++)
+    GmfSetLin(idx, GmfVertices, VerTab[i][0], VerTab[i][1], VerTab[i][2], dummyint);
+  
+  
+  /* ENFORCED VERTICES ========================== */
+  if (nbev) {
+    std::vector<std::vector<double> > ReqVerTab;
+    ReqVerTab.clear();
+    std::vector<double> aReqVerTab;
+    int solSize = 0;
+    for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
+      double x = vertexIt->first[0];
+      double y = vertexIt->first[1];
+      double z = vertexIt->first[2];
+      // Test if point is inside shape to mesh
+      gp_Pnt myPoint(x,y,z);
+      BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh());
+      scl.Perform(myPoint, 1e-7);
+      TopAbs_State result = scl.State();
+      if ( result != TopAbs_IN )
+        continue;
+      enfVertexSizes.push_back(vertexIt->second);
+      aReqVerTab.clear();
+      aReqVerTab.push_back(x);
+      aReqVerTab.push_back(y);
+      aReqVerTab.push_back(z);
+      ReqVerTab.push_back(aReqVerTab);
+      solSize++;
+    }
+
+    if (solSize) {
+      int idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+      if (!idxRequired)
+        return false;
+      int idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
+      if (!idxSol)
+        return false;
+      
+      int TypTab[] = {GmfSca};
+      GmfSetKwd(idxRequired, GmfVertices, solSize);
+      GmfSetKwd(idxSol, GmfSolAtVertices, solSize, 1, TypTab);
+      
+      for (int i=0;i<solSize;i++) {
+        double solTab[] = {enfVertexSizes.at(i)};
+        GmfSetLin(idxRequired, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
+        GmfSetLin(idxSol, GmfSolAtVertices, solTab);
+      }
+      GmfCloseMesh(idxRequired);
+      GmfCloseMesh(idxSol);
+    }
+  }
+
+  
+  /* ========================== FACES ========================== */
+  
+  int nbTriangles = 0, nbQuadrangles = 0, aSmdsID;
+  TopTools_IndexedMapOfShape facesMap, trianglesMap, quadranglesMap;
+  TIDSortedElemSet::const_iterator elemIt;
+  const SMESHDS_SubMesh* theSubMesh;
+  TopoDS_Shape aShape;
+  SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
+  const SMDS_MeshElement* aFace;
+  map<int,int>::const_iterator itOnMap;
+  std::vector<std::vector<int> > tt, qt,et;
+  tt.clear();
+  qt.clear();
+  et.clear();
+  std::vector<int> att, aqt, aet;
+  
+  TopExp::MapShapes( theMeshDS->ShapeToMesh(), TopAbs_FACE, facesMap );
+
+  for ( int i = 1; i <= facesMap.Extent(); ++i )
+    if (( theSubMesh  = theProxyMesh.GetSubMesh( facesMap(i))))
+    {
+      SMDS_ElemIteratorPtr it = theSubMesh->GetElements();
+      while (it->more())
+      {
+        const SMDS_MeshElement *elem = it->next();
+        if (elem->NbNodes() == 3)
+        {
+          trianglesMap.Add(facesMap(i));
+          nbTriangles ++;
+        }
+        else if (elem->NbNodes() == 4)
+        {
+          quadranglesMap.Add(facesMap(i));
+          nbQuadrangles ++;
+        }
+      }
+    }
+    
+  /* TRIANGLES ========================== */
+  if (nbTriangles) {
+    for ( int i = 1; i <= trianglesMap.Extent(); i++ )
+    {
+      aShape = trianglesMap(i);
+      theSubMesh = theProxyMesh.GetSubMesh(aShape);
+      if ( !theSubMesh ) continue;
+      itOnSubMesh = theSubMesh->GetElements();
+      while ( itOnSubMesh->more() )
+      {
+        aFace = itOnSubMesh->next();
+        itOnSubFace = aFace->nodesIterator();
+        att.clear();
+        while ( itOnSubFace->more() ) {
+          // find GHS3D ID
+          aSmdsID = itOnSubFace->next()->GetID();
+          itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
+          ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
+          att.push_back((*itOnMap).second);
+        }
+        tt.push_back(att);
+      }
+    }
+  }
+
+  if (theEnforcedTriangles.size()) {
+    // Iterate over the enforced triangles
+    for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
+      aFace = (*elemIt);
+      bool isOK = true;
+      itOnSubFace = aFace->nodesIterator();
+      att.clear();
+      while ( itOnSubFace->more() ) {
+        int aNodeID = itOnSubFace->next()->GetID();
+        itOnMap = theNodeId2NodeIndexMap.find(aNodeID);
+        if (itOnMap != theNodeId2NodeIndexMap.end())
+          att.push_back((*itOnMap).second);
+        else {
+          isOK = false;
+          theEnforcedTriangles.erase(elemIt);
+          break;
+        }
+      }
+      if (isOK)
+        tt.push_back(att);
+    }
+  }
+  
+  if (tt.size()) {
+    GmfSetKwd(idx, GmfTriangles, tt.size());
+    for (int i=0;i<tt.size();i++)
+      GmfSetLin(idx, GmfTriangles, tt[i][0], tt[i][1], tt[i][2], dummyint);
+  }
+    
+  /* QUADRANGLES ========================== */
+  if (nbQuadrangles) {
+    for ( int i = 1; i <= quadranglesMap.Extent(); i++ )
+    {
+      aShape = quadranglesMap(i);
+      theSubMesh = theProxyMesh.GetSubMesh(aShape);
+      if ( !theSubMesh ) continue;
+      itOnSubMesh = theSubMesh->GetElements();
+      while ( itOnSubMesh->more() )
+      {
+        aFace = itOnSubMesh->next();
+        itOnSubFace = aFace->nodesIterator();
+        aqt.clear();
+        while ( itOnSubFace->more() ) {
+          // find GHS3D ID
+          aSmdsID = itOnSubFace->next()->GetID();
+          itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
+          ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
+          aqt.push_back((*itOnMap).second);
+        }
+        qt.push_back(aqt);
+      }
+    }
+  }
+
+  if (theEnforcedQuadrangles.size()) {
+    // Iterate over the enforced triangles
+    for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) {
+      aFace = (*elemIt);
+      bool isOK = true;
+      itOnSubFace = aFace->nodesIterator();
+      aqt.clear();
+      while ( itOnSubFace->more() ) {
+        int aNodeID = itOnSubFace->next()->GetID();
+        itOnMap = theNodeId2NodeIndexMap.find(aNodeID);
+        if (itOnMap != theNodeId2NodeIndexMap.end())
+          aqt.push_back((*itOnMap).second);
+        else {
+          isOK = false;
+          theEnforcedQuadrangles.erase(elemIt);
+          break;
+        }
+      }
+      if (isOK)
+        qt.push_back(aqt);
+    }
+  }
+  
+  if (qt.size()) {
+    GmfSetKwd(idx, GmfQuadrilaterals, qt.size());
+    for (int i=0;i<qt.size();i++)
+      GmfSetLin(idx, GmfQuadrilaterals, qt[i][0], qt[i][1], qt[i][2], qt[i][3], dummyint);
+  }
+  
+
+  /* ========================== EDGES ========================== */
+
+  if (theEnforcedEdges.size()) {
+    // Iterate over the enforced edges
+    for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
+      aFace = (*elemIt);
+      bool isOK = true;
+      itOnSubFace = aFace->nodesIterator();
+      aet.clear();
+      while ( itOnSubFace->more() ) {
+        int aNodeID = itOnSubFace->next()->GetID();
+        itOnMap = theNodeId2NodeIndexMap.find(aNodeID);
+        if (itOnMap != theNodeId2NodeIndexMap.end())
+          aet.push_back((*itOnMap).second);
+        else {
+          isOK = false;
+          theEnforcedEdges.erase(elemIt);
+          break;
+        }
       }
+      if (isOK)
+        et.push_back(aet);
     }
   }
-  if ( !geomNormalOK)
-    return invalidID;
 
-  // compare normals
-  bool isReverse = ( meshNormal * geomNormal ) < 0;
-  if ( !isReverse )
-    return meshDS->ShapeToIndex( solid1 );
+  if (et.size()) {
+    GmfSetKwd(idx, GmfEdges, et.size());
+    for (int i=0;i<et.size();i++)
+      GmfSetLin(idx, GmfEdges, et[i][0], et[i][1], dummyint);
+  }
 
-  if ( solids.Extent() == 1 )
-    return HOLE_ID; // we are inside a hole
-  else
-    return meshDS->ShapeToIndex( solids(2) );
+  GmfCloseMesh(idx);
+  return true;
 }
 
+//=======================================================================
+//function : writeFaces
+//purpose  : Write Faces in case if generate 3D mesh with geometry
+//=======================================================================
+
+// static bool writeFaces (ofstream &             theFile,
+//                         const SMESH_ProxyMesh& theMesh,
+//                         const TopoDS_Shape&    theShape,
+//                         const map <int,int> &  theSmdsToGhs3dIdMap,
+//                         const map <int,int> &  theEnforcedNodeIdToGhs3dIdMap,
+//                         TIDSortedElemSet & theEnforcedEdges,
+//                         TIDSortedElemSet & theEnforcedTriangles,
+//                         TIDSortedElemSet & theEnforcedQuadrangles)
+// {
+//   // record structure:
+//   //
+//   // NB_ELEMS DUMMY_INT
+//   // Loop from 1 to NB_ELEMS
+//   // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
+// 
+//   TopoDS_Shape aShape;
+//   const SMESHDS_SubMesh* theSubMesh;
+//   const SMDS_MeshElement* aFace;
+//   const char* space    = "  ";
+//   const int   dummyint = 0;
+//   map<int,int>::const_iterator itOnMap;
+//   SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
+//   int nbNodes, aSmdsID;
+// 
+//   TIDSortedElemSet::const_iterator elemIt;
+//   int nbEnforcedEdges       = theEnforcedEdges.size();
+//   int nbEnforcedTriangles   = theEnforcedTriangles.size();
+//   int nbEnforcedQuadrangles = theEnforcedQuadrangles.size();
+//   // count triangles bound to geometry
+//   int nbTriangles = 0;
+// 
+//   TopTools_IndexedMapOfShape facesMap, trianglesMap, quadranglesMap;
+//   TopExp::MapShapes( theShape, TopAbs_FACE, facesMap );
+// 
+//   for ( int i = 1; i <= facesMap.Extent(); ++i )
+//     if (( theSubMesh  = theMesh.GetSubMesh( facesMap(i))))
+//       nbTriangles += theSubMesh->NbElements();
+// 
+//   std::cout << "    " << facesMap.Extent() << " shapes of 2D dimension and" << std::endl;
+//   if (nbEnforcedEdges+nbEnforcedTriangles+nbEnforcedQuadrangles)
+//     std::cout << "    " << nbEnforcedEdges+nbEnforcedTriangles+nbEnforcedQuadrangles 
+//                         << " enforced shapes:" << std::endl;
+//   if (nbEnforcedEdges)
+//     std::cout << "      " << nbEnforcedEdges << " enforced edges" << std::endl;
+//   if (nbEnforcedTriangles)
+//     std::cout << "      " << nbEnforcedTriangles << " enforced triangles" << std::endl;
+//   if (nbEnforcedQuadrangles)
+//     std::cout << "      " << nbEnforcedQuadrangles << " enforced quadrangles" << std::endl;
+//   std::cout << std::endl;
+// 
+// //   theFile << space << nbTriangles << space << dummyint << std::endl;
+//   std::ostringstream globalStream, localStream, aStream;
+// 
+//   //
+//   //        FACES : BEGIN
+//   //
+//   
+//   for ( int i = 1; i <= facesMap.Extent(); i++ )
+//   {
+//     aShape = facesMap(i);
+//     theSubMesh = theMesh.GetSubMesh(aShape);
+//     if ( !theSubMesh ) continue;
+//     itOnSubMesh = theSubMesh->GetElements();
+//     while ( itOnSubMesh->more() )
+//     {
+//       aFace = itOnSubMesh->next();
+//       nbNodes = aFace->NbNodes();
+// 
+//       localStream << nbNodes << space;
+// 
+//       itOnSubFace = aFace->nodesIterator();
+//       while ( itOnSubFace->more() ) {
+//         // find GHS3D ID
+//         aSmdsID = itOnSubFace->next()->GetID();
+//         itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
+//         // if ( itOnMap == theSmdsToGhs3dIdMap.end() ) {
+//         //   cout << "not found node: " << aSmdsID << endl;
+//         //   return false;
+//         // }
+//         ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
+// 
+//         localStream << (*itOnMap).second << space ;
+//       }
+// 
+//       // (NB_NODES + 1) times: DUMMY_INT
+//       for ( int j=0; j<=nbNodes; j++)
+//         localStream << dummyint << space ;
+// 
+//       localStream << std::endl;
+//     }
+//   }
+//   
+//   globalStream << localStream.str();
+//   localStream.str("");
+// 
+//   //
+//   //        FACES : END
+//   //
+// 
+//   //
+//   //        ENFORCED EDGES : BEGIN
+//   //
+//   
+//   // Iterate over the enforced edges
+//   int usedEnforcedEdges = 0;
+//   bool isOK;
+//   for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
+//     aFace = (*elemIt);
+//     isOK = true;
+//     itOnSubFace = aFace->nodesIterator();
+//     aStream.str("");
+//     aStream << "2" << space ;
+//     while ( itOnSubFace->more() ) {
+//       aSmdsID = itOnSubFace->next()->GetID();
+//       itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID);
+//       if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end())
+//         aStream << (*itOnMap).second << space;
+//       else {
+//         isOK = false;
+//         break;
+//       }
+//     }
+//     if (isOK) {
+//       for ( int j=0; j<=2; j++)
+//         aStream << dummyint << space ;
+// //       aStream << dummyint << space << dummyint;
+//       localStream << aStream.str() << std::endl;
+//       usedEnforcedEdges++;
+//     }
+//   }
+//   
+//   if (usedEnforcedEdges) {
+//     globalStream << localStream.str();
+//     localStream.str("");
+//   }
+// 
+//   //
+//   //        ENFORCED EDGES : END
+//   //
+//   //
+// 
+//   //
+//   //        ENFORCED TRIANGLES : BEGIN
+//   //
+//     // Iterate over the enforced triangles
+//   int usedEnforcedTriangles = 0;
+//   for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
+//     aFace = (*elemIt);
+//     isOK = true;
+//     itOnSubFace = aFace->nodesIterator();
+//     aStream.str("");
+//     aStream << "3" << space ;
+//     while ( itOnSubFace->more() ) {
+//       aSmdsID = itOnSubFace->next()->GetID();
+//       itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID);
+//       if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end())
+//         aStream << (*itOnMap).second << space;
+//       else {
+//         isOK = false;
+//         break;
+//       }
+//     }
+//     if (isOK) {
+//       for ( int j=0; j<=3; j++)
+//         aStream << dummyint << space ;
+//       localStream << aStream.str() << std::endl;
+//       usedEnforcedTriangles++;
+//     }
+//   }
+//   
+//   if (usedEnforcedTriangles) {
+//     globalStream << localStream.str();
+//     localStream.str("");
+//   }
+// 
+//   //
+//   //        ENFORCED TRIANGLES : END
+//   //
+// 
+//   //
+//   //        ENFORCED QUADRANGLES : BEGIN
+//   //
+//     // Iterate over the enforced quadrangles
+//   int usedEnforcedQuadrangles = 0;
+//   for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) {
+//     aFace = (*elemIt);
+//     isOK = true;
+//     itOnSubFace = aFace->nodesIterator();
+//     aStream.str("");
+//     aStream << "4" << space ;
+//     while ( itOnSubFace->more() ) {
+//       aSmdsID = itOnSubFace->next()->GetID();
+//       itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID);
+//       if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end())
+//         aStream << (*itOnMap).second << space;
+//       else {
+//         isOK = false;
+//         break;
+//       }
+//     }
+//     if (isOK) {
+//       for ( int j=0; j<=4; j++)
+//         aStream << dummyint << space ;
+//       localStream << aStream.str() << std::endl;
+//       usedEnforcedQuadrangles++;
+//     }
+//   }
+//   
+//   if (usedEnforcedQuadrangles) {
+//     globalStream << localStream.str();
+//     localStream.str("");
+//   }
+//   //
+//   //        ENFORCED QUADRANGLES : END
+//   //
+// 
+//   theFile
+//   << nbTriangles+usedEnforcedQuadrangles+usedEnforcedTriangles+usedEnforcedEdges
+//   << " 0" << std::endl
+//   << globalStream.str();
+// 
+//   return true;
+// }
+
+//=======================================================================
+//function : writePoints
+//purpose  : 
+//=======================================================================
+
+// static bool writePoints (ofstream &                       theFile,
+//                          SMESH_MesherHelper&              theHelper,
+//                          map <int,int> &                  theSmdsToGhs3dIdMap,
+//                          map <int,int> &                  theEnforcedNodeIdToGhs3dIdMap,
+//                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap,
+//                          GHS3DPlugin_Hypothesis::TID2SizeMap & theNodeIDToSizeMap,
+//                          GHS3DPlugin_Hypothesis::TEnforcedVertexValues & theEnforcedVertices,
+//                          TIDSortedNodeSet & theEnforcedNodes)
+// {
+//   // record structure:
+//   //
+//   // NB_NODES
+//   // Loop from 1 to NB_NODES
+//   //   X Y Z DUMMY_INT
+// 
+//   SMESHDS_Mesh * theMesh = theHelper.GetMeshDS();
+//   int nbNodes = theMesh->NbNodes();
+//   if ( nbNodes == 0 )
+//     return false;
+//   int nbEnforcedVertices = theEnforcedVertices.size();
+//   int nbEnforcedNodes    = theEnforcedNodes.size();
+// 
+//   int aGhs3dID = 1;
+//   SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
+//   const SMDS_MeshNode* node;
+// 
+//   // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
+//   // The problem is in nodes on degenerated edges, we need to skip them
+//   if ( theHelper.HasDegeneratedEdges() )
+//   {
+//     // here we decrease total nb of nodes by nb of nodes on degenerated edges
+//     set<int> checkedSM;
+//     for (TopExp_Explorer e(theMesh->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
+//     {
+//       SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
+//       if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() )) {
+//         if ( sm->GetSubMeshDS() )
+//           nbNodes -= sm->GetSubMeshDS()->NbNodes();
+//       }
+//     }
+//   }
+// 
+//   const bool isQuadMesh = 
+//     theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
+//     theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
+//     theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC );
+//   if ( isQuadMesh )
+//   {
+//     // descrease nbNodes by nb of medium nodes
+//     while ( it->more() )
+//     {
+//       node = it->next();
+//       if ( !theHelper.IsDegenShape( node->getshapeId() ))
+//         nbNodes -= int( theHelper.IsMedium( node ));
+//     }
+//     it = theMesh->nodesIterator();
+//   }
+// 
+//   const char* space    = "  ";
+//   const int   dummyint = 0;
+// 
+//   // NB_NODES
+//   std::cout << std::endl;
+//   std::cout << "The initial 2D mesh contains :" << std::endl;
+//   std::cout << "    " << nbNodes << " nodes" << std::endl;
+//   if (nbEnforcedVertices > 0)
+//     std::cout << "    " << nbEnforcedVertices << " enforced vertices" << std::endl;
+//   if (nbEnforcedNodes > 0)
+//     std::cout << "    " << nbEnforcedNodes << " enforced nodes" << std::endl;
+// 
+// //   std::cout << std::endl;
+// //   std::cout << "Start writing in 'points' file ..." << std::endl;
+//   theFile << nbNodes << space << std::endl;
+// 
+//   // Loop from 1 to NB_NODES
+// 
+//   while ( it->more() )
+//   {
+//     node = it->next();
+//     if (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238
+//         theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674
+//       continue;
+// 
+//     theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
+//     theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
+//     aGhs3dID++;
+// 
+//     // X Y Z DUMMY_INT
+//     theFile
+//     << node->X() << space 
+//     << node->Y() << space 
+//     << node->Z() << space 
+//     << dummyint << space ;
+//     theFile << std::endl;
+// 
+//   }
+//   
+//   // Iterate over the enforced nodes
+//   TIDSortedNodeSet::const_iterator nodeIt;
+//   std::map<int,double> enfVertexIndexSizeMap;
+//   if (nbEnforcedNodes) {
+//     for(nodeIt = theEnforcedNodes.begin() ; nodeIt != theEnforcedNodes.end() ; ++nodeIt) {
+//       double x = (*nodeIt)->X();
+//       double y = (*nodeIt)->Y();
+//       double z = (*nodeIt)->Z();
+//       // Test if point is inside shape to mesh
+//       gp_Pnt myPoint(x,y,z);
+//       BRepClass3d_SolidClassifier scl(theMesh->ShapeToMesh());
+//       scl.Perform(myPoint, 1e-7);
+//       TopAbs_State result = scl.State();
+//       if ( result != TopAbs_IN )
+//         continue;
+//       std::vector<double> coords;
+//       coords.push_back(x);
+//       coords.push_back(y);
+//       coords.push_back(z);
+//       if (theEnforcedVertices.find(coords) != theEnforcedVertices.end())
+//         continue;
+//         
+//       double size = theNodeIDToSizeMap.find((*nodeIt)->GetID())->second;
+//   //       theGhs3dIdToNodeMap.insert( make_pair( nbNodes + i, (*nodeIt) ));
+//   //       MESSAGE("Adding enforced node (" << x << "," << y <<"," << z << ")");
+//       // X Y Z PHY_SIZE DUMMY_INT
+//       theFile
+//       << x << space 
+//       << y << space 
+//       << z << space
+//       << size << space
+//       << dummyint << space;
+//       theFile << std::endl;
+//       theEnforcedNodeIdToGhs3dIdMap.insert( make_pair( (*nodeIt)->GetID(), aGhs3dID ));
+//       enfVertexIndexSizeMap[aGhs3dID] = -1;
+//       aGhs3dID++;
+//   //     else
+//   //         MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added ");
+//     }
+//   }
+//   
+//   if (nbEnforcedVertices) {
+//     // Iterate over the enforced vertices
+//     GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
+// //     int i = 1;
+//     for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
+//       double x = vertexIt->first[0];
+//       double y = vertexIt->first[1];
+//       double z = vertexIt->first[2];
+//       // Test if point is inside shape to mesh
+//       gp_Pnt myPoint(x,y,z);
+//       BRepClass3d_SolidClassifier scl(theMesh->ShapeToMesh());
+//       scl.Perform(myPoint, 1e-7);
+//       TopAbs_State result = scl.State();
+//       if ( result != TopAbs_IN )
+//         continue;
+//   //         MESSAGE("Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second);
+//           // X Y Z PHY_SIZE DUMMY_INT
+//       theFile
+//       << x << space 
+//       << y << space 
+//       << z << space
+//       << vertexIt->second << space 
+//       << dummyint << space;
+//       theFile << std::endl;
+//       enfVertexIndexSizeMap[aGhs3dID] = vertexIt->second;
+//       aGhs3dID++;
+//       
+//   //     else
+//   //         MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added ");
+//     }
+//   }
+//   theFile << std::endl;
+//   
+//   
+// //   std::cout << std::endl;
+// //   std::cout << "End writing in 'points' file." << std::endl;
+// 
+//   return true;
+// }
+
+
 //=======================================================================
 //function : readResultFile
 //purpose  : 
@@ -830,12 +2446,17 @@ static bool readResultFile(const int                       fileOpen,
                            GHS3DPlugin_GHS3D*              theAlgo,
 #endif
                            SMESH_MesherHelper&             theHelper,
+//                            SMESH_Mesh&                     theMesh,
                            TopoDS_Shape                    tabShape[],
                            double**                        tabBox,
                            const int                       nbShape,
                            map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
                            bool                            toMeshHoles,
-                           int                             nbEnforcedVertices)
+                           int                             nbEnforcedVertices,
+                           int                             nbEnforcedNodes,
+                           TIDSortedElemSet &              theEnforcedEdges,
+                           TIDSortedElemSet &              theEnforcedTriangles,
+                           TIDSortedElemSet &              theEnforcedQuadrangles)
 {
   MESSAGE("GHS3DPlugin_GHS3D::readResultFile()");
   Kernel_Utils::Localizer loc;
@@ -917,7 +2538,7 @@ static bool readResultFile(const int                       fileOpen,
     for (int iCoor=0; iCoor < 3; iCoor++)
       coord[ iCoor ] = strtod(ptr, &ptr);
     nodeAssigne[ iNode ] = 1;
-    if ( iNode > (nbInputNodes-nbEnforcedVertices) ) {
+    if ( iNode > (nbInputNodes-(nbEnforcedVertices+nbEnforcedNodes)) ) {
       // Creating SMESH nodes
       // - for enforced vertices
       // - for vertices of forced edges
@@ -950,6 +2571,7 @@ static bool readResultFile(const int                       fileOpen,
       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
       try {
         OCC_CATCH_SIGNALS;
+//         tabID[i] = findShapeID( theHelper, n1, n2, n3, toMeshHoles );
         tabID[i] = findShapeID( *theHelper.GetMesh(), n1, n2, n3, toMeshHoles );
         // -- 0020330: Pb with ghs3d as a submesh
         // check that found shape is to be meshed
@@ -1094,174 +2716,10 @@ static bool readResultFile(const int                       fileOpen,
   return true;
 }
 
-//=======================================================================
-//function : readResultFile
-//purpose  : 
-//=======================================================================
-
-static bool readResultFile(const int                      fileOpen,
-#ifdef WNT
-                           const char*                    fileName,
-#endif
-#ifdef WITH_SMESH_CANCEL_COMPUTE
-                           GHS3DPlugin_GHS3D*             theAlgo,
-#endif
-                           SMESH_MesherHelper&            theHelper,
-                           TopoDS_Shape                   aSolid,
-                           vector <const SMDS_MeshNode*>& theNodeByGhs3dId,
-                           int                            nbEnforcedVertices)
-{
-  SMESHDS_Mesh* theMeshDS = theHelper.GetMeshDS();
-
-  Kernel_Utils::Localizer loc;
-  struct stat  status;
-  size_t       length;
-
-  char *ptr, *mapPtr;
-  char *tetraPtr;
-  char *shapePtr;
-
-  int fileStat;
-  int nbElems, nbNodes, nbInputNodes;
-  int nodeId, triangleId;
-  int nbTriangle;
-  int ID, shapeID;
-
-  int *tab;
-  double *coord;
-  const SMDS_MeshNode **node;
-
-  tab   = new int[3];
-  coord = new double[3];
-  node  = new const SMDS_MeshNode*[4];
-
-  SMDS_MeshNode * aNewNode;
-  map <int,const SMDS_MeshNode*>::iterator IdNode;
-  SMDS_MeshElement* aTet;
-
-  // Read the file state
-  fileStat = fstat(fileOpen, &status);
-  length   = status.st_size;
-
-  // Mapping the result file into memory
-#ifdef WNT
-  HANDLE fd = CreateFile(fileName, GENERIC_READ, FILE_SHARE_READ,
-                         NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
-  HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
-                                        0, (DWORD)length, NULL);
-  ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
-#else
-  ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
-#endif
-  mapPtr = ptr;
-
-  ptr      = readMapIntLine(ptr, tab);
-  tetraPtr = ptr;
-
-  nbElems      = tab[0];
-  nbNodes      = tab[1];
-  nbInputNodes = tab[2];
-
-  theNodeByGhs3dId.resize( nbNodes );
-
-  // Reading the nodeId
-  for (int i=0; i < 4*nbElems; i++)
-    nodeId = strtol(ptr, &ptr, 10);
-
-  // Issue 0020682. Avoid creating nodes and tetras at place where
-  // volumic elements already exist
-  SMESH_ElementSearcher* elemSearcher = 0;
-  vector< const SMDS_MeshElement* > foundVolumes;
-  if ( theHelper.GetMesh()->NbVolumes() > 0 )
-    elemSearcher = SMESH_MeshEditor( theHelper.GetMesh() ).GetElementSearcher();
-
-  // Reading the nodeCoord and update the nodeMap
-  shapeID = theMeshDS->ShapeToIndex( aSolid );
-  for (int iNode=0; iNode < nbNodes; iNode++) {
-#ifdef WITH_SMESH_CANCEL_COMPUTE
-    if(theAlgo->computeCanceled())
-      return false;
-#endif
-    for (int iCoor=0; iCoor < 3; iCoor++)
-      coord[ iCoor ] = strtod(ptr, &ptr);
-    if ((iNode+1) > (nbInputNodes-nbEnforcedVertices)) {
-      // Issue 0020682. Avoid creating nodes and tetras at place where
-      // volumic elements already exist
-      if ( elemSearcher &&
-           elemSearcher->FindElementsByPoint( gp_Pnt(coord[0],coord[1],coord[2]),
-                                              SMDSAbs_Volume, foundVolumes ))
-      {
-        theNodeByGhs3dId[ iNode ] = 0;
-      }
-      else
-      {
-        aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
-        theMeshDS->SetNodeInVolume( aNewNode, shapeID );
-        theNodeByGhs3dId[ iNode ] = aNewNode;
-      }
-    }
-  }
-
-  // Reading the triangles
-  nbTriangle = strtol(ptr, &ptr, 10);
-
-  for (int i=0; i < 3*nbTriangle; i++)
-    triangleId = strtol(ptr, &ptr, 10);
-
-  shapePtr = ptr;
-
-  if ( theHelper.IsQuadraticMesh() != SMESH_MesherHelper::LINEAR )
-    theHelper.IsQuadraticSubMesh( aSolid );
-
-  // Associating the tetrahedrons to the shapes
-
-  for (int iElem = 0; iElem < nbElems; iElem++) {
-#ifdef WITH_SMESH_CANCEL_COMPUTE
-    if(theAlgo->computeCanceled())
-      return false;
-#endif
-    for (int iNode = 0; iNode < 4; iNode++) {
-      ID = strtol(tetraPtr, &tetraPtr, 10);
-      node[ iNode ] = theNodeByGhs3dId[ ID-1 ];
-    }
-    if ( elemSearcher )
-    {
-      // Issue 0020682. Avoid creating nodes and tetras at place where
-      // volumic elements already exist
-      if ( !node[1] || !node[0] || !node[2] || !node[3] )
-        continue;
-      if ( elemSearcher->FindElementsByPoint(( SMESH_TNodeXYZ(node[0]) +
-                                               SMESH_TNodeXYZ(node[1]) +
-                                               SMESH_TNodeXYZ(node[2]) +
-                                               SMESH_TNodeXYZ(node[3]) ) / 4.,
-                                             SMDSAbs_Volume, foundVolumes ))
-        continue;
-    }
-    aTet = theHelper.AddVolume( node[1], node[0], node[2], node[3],
-                                /*id=*/0, /*force3d=*/true);
-    theMeshDS->SetMeshElementOnShape( aTet, shapeID );
-  }
-  if ( nbElems )
-    cout << nbElems << " tetrahedrons have been associated to " << nbTriangle << " shapes" << endl;
-#ifdef WNT
-  UnmapViewOfFile(mapPtr);
-  CloseHandle(hMapObject);
-  CloseHandle(fd);
-#else
-  munmap(mapPtr, length);
-#endif
-  close(fileOpen);
-
-  delete [] tab;
-  delete [] coord;
-  delete [] node;
-
-  return true;
-}
 
 //=============================================================================
 /*!
- *Here we are going to use the GHS3D mesher
+ *Here we are going to use the GHS3D mesher with geometry
  */
 //=============================================================================
 
@@ -1305,39 +2763,66 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
   TCollection_AsciiString aGenericName
     = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str();
 
-  TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
-  TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
-  aFacesFileName  = aGenericName + ".faces";  // in faces
-  aPointsFileName = aGenericName + ".points"; // in points
-  aResultFileName = aGenericName + ".noboite";// out points and volumes
-  aBadResFileName = aGenericName + ".boite";  // out bad result
-  aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
-  aLogFileName    = aGenericName + ".log";    // log
-
-  // -----------------
-  // make input files
-  // -----------------
-
-  ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
-  ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
+  TCollection_AsciiString aResultFileName;
+  TCollection_AsciiString aLogFileName    = aGenericName + ".log";    // log
+// #if GHS3D_VERSION < 42
+//   TCollection_AsciiString aFacesFileName, aPointsFileName;
+//   TCollection_AsciiString aBadResFileName, aBbResFileName;
+//   aFacesFileName  = aGenericName + ".faces";  // in faces
+//   aPointsFileName = aGenericName + ".points"; // in points
+//   aResultFileName = aGenericName + ".noboite";// out points and volumes
+//   aBadResFileName = aGenericName + ".boite";  // out bad result
+//   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
+// 
+//   // -----------------
+//   // make input files
+//   // -----------------
+// 
+//   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
+//   ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
+// 
+//   Ok =
+//     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
+//   if (!Ok) {
+//     INFOS( "Can't write into " << aFacesFileName);
+//     return error(SMESH_Comment("Can't write into ") << aFacesFileName);
+//   }
+// #else
+  TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName;
+  TCollection_AsciiString aResultGMFFileName;
 
-  Ok =
-    aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
-  if (!Ok) {
-    INFOS( "Can't write into " << aFacesFileName);
-    return error(SMESH_Comment("Can't write into ") << aFacesFileName);
-  }
-  map <int,int> aSmdsToGhs3dIdMap;
+#ifdef _DEBUG_
+  aGMFFileName    = aGenericName + ".mesh"; // GMF mesh file
+  // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
+  aResultGMFFileName = aGenericName + "Vol.mesh"; // GMF mesh file
+  aResultFileName = aGenericName + ".noboite";// out points and volumes
+//   aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file
+  aRequiredVerticesFileName    = aGenericName + "_required.mesh"; // GMF required vertices mesh file
+  aSolFileName    = aGenericName + "_required.sol"; // GMF solution file
+#else
+  aGMFFileName    = aGenericName + ".mesh"; // GMF mesh file
+//   aGMFFileName    = aGenericName + ".meshb"; // GMF mesh file
+  // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
+  aResultGMFFileName = aGenericName + "Vol.meshb"; // GMF mesh file
+  aResultFileName = aGenericName + ".noboite";// out points and volumes
+//   aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file
+  aRequiredVerticesFileName    = aGenericName + "_required.meshb"; // GMF required vertices mesh file
+  aSolFileName    = aGenericName + "_required.solb"; // GMF solution file
+#endif
+  map <int,int> aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap;
   map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
-  map<vector<double>,double> enforcedVertices;
-  int nbEnforcedVertices = 0;
-  try {
-    enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
-    nbEnforcedVertices = enforcedVertices.size();
-  }
-  catch(...) {
-  }
-
+  std::map <int, int> nodeID2nodeIndexMap;
+  GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
+  TIDSortedNodeSet enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp);
+  TIDSortedElemSet enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp);
+  TIDSortedElemSet enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp);
+  TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp);
+  GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp);
+//   GHS3DPlugin_Hypothesis::TID2SizeMap elemIDToSizeMap = GHS3DPlugin_Hypothesis::GetElementIDToSizeMap(_hyp);
+
+  int nbEnforcedVertices = enforcedVertices.size();
+  int nbEnforcedNodes = enforcedNodes.size();
+  
   SMESH_MesherHelper helper( theMesh );
   helper.SetSubShape( theShape );
 
@@ -1369,18 +2854,29 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
       if ( !proxyMesh )
         return false;
     }
-
-    Ok = (writePoints( aPointsFile, helper, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, enforcedVertices)
-          &&
-          writeFaces ( aFacesFile, *proxyMesh, theShape, aSmdsToGhs3dIdMap ));
+// #if GHS3D_VERSION < 42
+//     Ok = (writePoints( aPointsFile, helper, 
+//                        aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, aGhs3dIdToNodeMap, 
+//                        nodeIDToSizeMap,
+//                        enforcedVertices, enforcedNodes)
+//           &&
+//           writeFaces ( aFacesFile, *proxyMesh, theShape, 
+//                        aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap,
+//                        enforcedEdges, enforcedTriangles, enforcedQuadrangles));
+// #else
+    Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(),
+                      helper, *proxyMesh,
+                      aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap,
+                      enforcedNodes, enforcedEdges, enforcedTriangles, enforcedQuadrangles,
+                      enforcedVertices);
+// #endif
   }
 
   // Write aSmdsToGhs3dIdMap to temp file
   TCollection_AsciiString aSmdsToGhs3dIdMapFileName;
   aSmdsToGhs3dIdMapFileName = aGenericName + ".ids";  // ids relation
   ofstream aIdsFile  ( aSmdsToGhs3dIdMapFileName.ToCString()  , ios::out);
-  Ok =
-    aIdsFile.rdbuf()->is_open();
+  Ok = aIdsFile.rdbuf()->is_open();
   if (!Ok) {
     INFOS( "Can't write into " << aSmdsToGhs3dIdMapFileName);
     return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName);
@@ -1391,14 +2887,22 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
     aIdsFile << myit->first << " " << myit->second << std::endl;
   }
 
-  aFacesFile.close();
-  aPointsFile.close();
   aIdsFile.close();
+// #if GHS3D_VERSION < 42
+//   aFacesFile.close();
+//   aPointsFile.close();
+// #endif
   
   if ( ! Ok ) {
     if ( !_keepFiles ) {
-      removeFile( aFacesFileName );
-      removeFile( aPointsFileName );
+// #if GHS3D_VERSION < 42
+//       removeFile( aFacesFileName );
+//       removeFile( aPointsFileName );
+// #else
+      removeFile( aGMFFileName );
+      removeFile( aRequiredVerticesFileName );
+      removeFile( aSolFileName );
+// #endif
       removeFile( aSmdsToGhs3dIdMapFileName );
     }
     return error(COMPERR_BAD_INPUT_MESH);
@@ -1409,9 +2913,14 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
   // run ghs3d mesher
   // -----------------
 
-  TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
+  TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
+  // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
-  cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
+  cmd += TCollection_AsciiString(" -IM ");
+//   cmd += TCollection_AsciiString(" --in ") + aGenericName;
+//   cmd += TCollection_AsciiString(" --required_vertices ") + aRequiredVerticesFileName;
+//    cmd += TCollection_AsciiString(" --out ") + aGenericName;
+  cmd += TCollection_AsciiString(" -Om 1>" ) + aLogFileName;  // dump into file
 
   std::cout << std::endl;
   std::cout << "Ghs3d execution..." << std::endl;
@@ -1430,6 +2939,7 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
   // read a result
   // --------------
 
+// #if GHS3D_VERSION < 42
   // Mapping the result file
 
   int fileOpen;
@@ -1454,9 +2964,37 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 #ifdef WITH_SMESH_CANCEL_COMPUTE
                          this,
 #endif
-                         helper, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
-                         toMeshHoles, nbEnforcedVertices );
+                         /*theMesh, */helper, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
+                         toMeshHoles, 
+                         nbEnforcedVertices, nbEnforcedNodes, 
+                         enforcedEdges, enforcedTriangles, enforcedQuadrangles );
   }
+// /*/*#else
+// #ifndef GMF_HAS_SUBDOMAIN_INFO
+//   // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
+//   
+//   int fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
+//   if ( fileOpen < 0 ) {
+//     std::cout << std::endl;
+//     std::cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << std::endl;
+//     std::cout << "Log: " << aLogFileName << std::endl;
+//     Ok = false;
+//   }
+//   else {
+// #endif
+//       Ok = readGMFFile(
+// #ifndef GMF_HAS_SUBDOMAIN_INFO
+//                        fileOpen,
+// #endif
+//                        aGenericName.ToCString(), theMesh,
+//                        _nbShape, tabShape, tabBox, 
+//                        aGhs3dIdToNodeMap, toMeshHoles,
+//                        nbEnforcedVertices, nbEnforcedNodes, 
+//                        enforcedNodes, enforcedTriangles, enforcedQuadrangles);
+// #ifndef GMF_HAS_SUBDOMAIN_INFO
+//   }
+// #endif
+// #endif*/*/
 
   // ---------------------
   // remove working files
@@ -1482,17 +3020,21 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
   }
 
   if ( !_keepFiles ) {
+// #if GHS3D_VERSION < 42
+//     removeFile( aFacesFileName );
+//     removeFile( aPointsFileName );
+//     removeFile( aResultFileName );
+//     removeFile( aBadResFileName );
+//     removeFile( aBbResFileName );
+// #endif
+    removeFile( aSmdsToGhs3dIdMapFileName );
+  // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2)
+
 #ifdef WITH_SMESH_CANCEL_COMPUTE
     if (! Ok)
       if(_compute_canceled)
         removeFile( aLogFileName );
 #endif
-    removeFile( aFacesFileName );
-    removeFile( aPointsFileName );
-    removeFile( aResultFileName );
-    removeFile( aBadResFileName );
-    removeFile( aBbResFileName );
-    removeFile( aSmdsToGhs3dIdMapFileName );
   }
   std::cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
   if ( !Ok )
@@ -1525,37 +3067,53 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
   TCollection_AsciiString aGenericName
     = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str();
 
-  TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
-  TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
-  aFacesFileName  = aGenericName + ".faces";  // in faces
-  aPointsFileName = aGenericName + ".points"; // in points
-  aResultFileName = aGenericName + ".noboite";// out points and volumes
-  aBadResFileName = aGenericName + ".boite";  // out bad result
-  aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
-  aLogFileName    = aGenericName + ".log";    // log
-
-  // -----------------
-  // make input files
-  // -----------------
-
-  ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
-  ofstream aPointsFile  ( aPointsFileName.ToCString()  , ios::out);
-  bool Ok =
-    aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
-
-  if (!Ok)
-    return error( SMESH_Comment("Can't write into ") << aPointsFileName);
+  TCollection_AsciiString aLogFileName    = aGenericName + ".log";    // log
+  TCollection_AsciiString aResultFileName;
+  bool Ok;
+// #if GHS3D_VERSION < 42
+//   TCollection_AsciiString aFacesFileName, aPointsFileName;
+//   TCollection_AsciiString aBadResFileName, aBbResFileName;
+//   aFacesFileName  = aGenericName + ".faces";  // in faces
+//   aPointsFileName = aGenericName + ".points"; // in points
+//   aResultFileName = aGenericName + ".noboite";// out points and volumes
+//   aBadResFileName = aGenericName + ".boite";  // out bad result
+//   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
+// 
+//   // -----------------
+//   // make input files
+//   // -----------------
+// 
+//   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
+//   ofstream aPointsFile  ( aPointsFileName.ToCString()  , ios::out);
+//   Ok = aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
+//   if (!Ok) {
+//     INFOS( "Can't write into " << aFacesFileName);
+//     return error( SMESH_Comment("Can't write into ") << aFacesFileName);
+//   }
+// #else
+  TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName;
+#ifdef _DEBUG_
+  aGMFFileName    = aGenericName + ".mesh"; // GMF mesh file
+  aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file
+  aRequiredVerticesFileName    = aGenericName + "_required.mesh"; // GMF required vertices mesh file
+  aSolFileName    = aGenericName + "_required.sol"; // GMF solution file
+#else
+  aGMFFileName    = aGenericName + ".meshb"; // GMF mesh file
+  aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file
+  aRequiredVerticesFileName    = aGenericName + "_required.meshb"; // GMF required vertices mesh file
+  aSolFileName    = aGenericName + ".solb"; // GMF solution file
+#endif
+// #endif
   
-  GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices;
-  int nbEnforcedVertices = 0;
-  try {
-    enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
-    nbEnforcedVertices = enforcedVertices.size();
-  }
-  catch(...) {
-  }
-
-  vector <const SMDS_MeshNode*> aNodeByGhs3dId;
+  std::map <int, int> nodeID2nodeIndexMap;
+  GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
+  TIDSortedNodeSet enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp);
+  TIDSortedElemSet enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp);
+  TIDSortedElemSet enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp);
+  TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp);
+  GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp);
+
+  vector <const SMDS_MeshNode*> aNodeByGhs3dId, anEnforcedNodeByGhs3dId;
   {
     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
     if ( theMesh.NbQuadrangles() > 0 )
@@ -1564,30 +3122,57 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
       aQuad2Trias->Compute( theMesh );
       proxyMesh.reset( aQuad2Trias );
     }
-
-    Ok = (writeFaces ( aFacesFile, *proxyMesh, aNodeByGhs3dId ) &&
-          writePoints( aPointsFile, &theMesh, aNodeByGhs3dId,enforcedVertices));
-  }  
-  aFacesFile.close();
-  aPointsFile.close();
-  
-  if ( ! Ok ) {
-    if ( !_keepFiles ) {
-      removeFile( aFacesFileName );
-      removeFile( aPointsFileName );
-    }
-    return error(COMPERR_BAD_INPUT_MESH);
+// #if GHS3D_VERSION < 42
+//     Ok = (writeFaces ( aFacesFile, *proxyMesh, &theMesh, aNodeByGhs3dId, anEnforcedNodeByGhs3dId,
+//                        enforcedEdges, enforcedTriangles, enforcedQuadrangles ) &&
+//           writePoints( aPointsFile, &theMesh, aNodeByGhs3dId, anEnforcedNodeByGhs3dId,
+//                        nodeIDToSizeMap, enforcedVertices, enforcedNodes));
+//     int nbEnforcedVertices = enforcedVertices.size();
+//     int nbEnforcedNodes = enforcedNodes.size();
+// #else
+    Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(),
+                      *proxyMesh, &theMesh,
+                      aNodeByGhs3dId, anEnforcedNodeByGhs3dId,
+                      enforcedNodes, enforcedEdges, enforcedTriangles, enforcedQuadrangles,
+                      enforcedVertices);
+// #endif
   }
-  removeFile( aResultFileName ); // needed for boundary recovery module usage
+  
+  TIDSortedNodeSet enforcedNodesFromEnforcedElem;
+  for (int i=0;i<anEnforcedNodeByGhs3dId.size();i++)
+    enforcedNodesFromEnforcedElem.insert(anEnforcedNodeByGhs3dId[i]);
+
+// #if GHS3D_VERSION < 42
+//   aFacesFile.close();
+//   aPointsFile.close();
+//   
+//   if ( ! Ok ) {
+//     if ( !_keepFiles ) {
+//       removeFile( aFacesFileName );
+//       removeFile( aPointsFileName );
+//     }
+//     return error(COMPERR_BAD_INPUT_MESH);
+//   }
+//   removeFile( aResultFileName ); // needed for boundary recovery module usage
+// #endif
 
   // -----------------
   // run ghs3d mesher
   // -----------------
 
-  TCollection_AsciiString cmd =
-    (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str();
-  cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
+  TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str());
+// #if GHS3D_VERSION < 42
+//   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
+// #else
+  cmd += TCollection_AsciiString(" --in ") + aGenericName;
+//   cmd += TCollection_AsciiString(" --required_vertices ") + aRequiredVerticesFileName;
+  cmd += TCollection_AsciiString(" --out ") + aResultFileName;
+// #endif
   cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
+
+  std::cout << std::endl;
+  std::cout << "Ghs3d execution..." << std::endl;
+  std::cout << cmd << std::endl;
   
 #ifdef WITH_SMESH_CANCEL_COMPUTE
   _compute_canceled = false;
@@ -1595,28 +3180,41 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 
   system( cmd.ToCString() ); // run
 
+  std::cout << std::endl;
+  std::cout << "End of Ghs3d execution !" << std::endl;
+
   // --------------
   // read a result
   // --------------
-  int fileOpen;
-  fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
-  if ( fileOpen < 0 ) {
-    std::cout << std::endl;
-    std::cout << "Error when opening the " << aResultFileName.ToCString() << " file" << std::endl;
-    std::cout << "Log: " << aLogFileName << std::endl;
-    std::cout << std::endl;
-    Ok = false;
-  }
-  else {
-    Ok = readResultFile( fileOpen,
-#ifdef WNT
-                         aResultFileName.ToCString(),
-#endif
+// #if GHS3D_VERSION < 42
+//   int fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
+//   if ( fileOpen < 0 ) {
+//     std::cout << std::endl;
+//     std::cout << "Error when opening the " << aResultFileName.ToCString() << " file" << std::endl;
+//     std::cout << "Log: " << aLogFileName << std::endl;
+//     std::cout << std::endl;
+//     Ok = false;
+//   }
+//   else {
+//     Ok = readResultFile( fileOpen,
+// #ifdef WNT
+//                          aResultFileName.ToCString(),
+// #endif
+// #ifdef WITH_SMESH_CANCEL_COMPUTE
+//                          this,
+// #endif
+//                          theMesh, theShape ,aNodeByGhs3dId, anEnforcedNodeByGhs3dId,
+//                          nbEnforcedVertices, nbEnforcedNodes, 
+//                          enforcedEdges, enforcedTriangles, enforcedQuadrangles );
+//   }
+// #else
+  Ok = readGMFFile(aResultFileName.ToCString(),
 #ifdef WITH_SMESH_CANCEL_COMPUTE
-                         this,
+                   this,
 #endif
-                         *theHelper, theShape ,aNodeByGhs3dId, nbEnforcedVertices );
-  }
+                   theHelper, 
+                   enforcedNodesFromEnforcedElem, enforcedTriangles, enforcedQuadrangles);
+// #endif
   
   // ---------------------
   // remove working files
@@ -1639,7 +3237,7 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
   }
-
+// #if GHS3D_VERSION < 42
   if ( !_keepFiles )
   {
 #ifdef WITH_SMESH_CANCEL_COMPUTE
@@ -1647,13 +3245,13 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
       if(_compute_canceled)
         removeFile( aLogFileName );
 #endif
-    removeFile( aFacesFileName );
-    removeFile( aPointsFileName );
-    removeFile( aResultFileName );
-    removeFile( aBadResFileName );
-    removeFile( aBbResFileName );
+//     removeFile( aFacesFileName );
+//     removeFile( aPointsFileName );
+//     removeFile( aResultFileName );
+//     removeFile( aBadResFileName );
+//     removeFile( aBbResFileName );
   }
-  
+// #endif
   return Ok;
 }
 
@@ -2253,3 +3851,14 @@ bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh,
   return true;
 }
 
+bool GHS3DPlugin_GHS3D::importGMFMesh(const char* theGMFFileName, SMESH_Mesh& theMesh)
+{
+  SMESH_MesherHelper* helper = new SMESH_MesherHelper(theMesh );
+  TIDSortedElemSet dummyElemSet;
+  TIDSortedNodeSet dummyNodeSet;
+  return readGMFFile(theGMFFileName, 
+#ifdef WITH_SMESH_CANCEL_COMPUTE
+                   this,
+#endif
+                   helper, dummyNodeSet , dummyElemSet, dummyElemSet);
+}
index 1c2c32082a5c0272565681e99be7ad56f2bc5c89..3888942eb37be55059d1483a20300ff7ed93db6e 100644 (file)
 #define _GHS3DPlugin_GHS3D_HXX_
 
 #include "SMESH_3D_Algo.hxx"
+#include "SMESH_Gen.hxx"
 
 #include <map>
 #include <vector>
 
+extern "C"
+{
+  #include "libmesh5.h"
+}
+
+#ifndef GMFVERSION
+#define GMFVERSION GmfDouble
+#endif
+#define GMFDIMENSION 3
+
 class GHS3DPlugin_Hypothesis;
 class SMDS_MeshNode;
 class SMESH_Mesh;
@@ -62,6 +73,8 @@ public:
   virtual bool Compute(SMESH_Mesh&         theMesh,
                        SMESH_MesherHelper* aHelper);
 
+  bool importGMFMesh(const char* aGMFFileName, SMESH_Mesh& aMesh);
+
 private:
 
   bool storeErrorDescription(const TCollection_AsciiString& logFile,
index c63525a0cd92e5a9ed4810172710cf04703d56a2..cee96775f74b40a058a7393391a711050edc6225 100644 (file)
 //
 #include "GHS3DPlugin_GHS3D_i.hxx"
 #include "SMESH_Gen.hxx"
+#include "SMESH_Mesh_i.hxx"
+#include "SMESH_Gen_i.hxx"
 #include "GHS3DPlugin_GHS3D.hxx"
+#include "SMESH_PythonDump.hxx"
 
 #include "utilities.h"
+#include <cstring>
 
 using namespace std;
 
@@ -80,3 +84,36 @@ GHS3DPlugin_GHS3D_i::~GHS3DPlugin_GHS3D_i()
   return ( ::GHS3DPlugin_GHS3D* )myBaseImpl;
 }
 
+//=============================================================================
+/*!
+ *  GHS3DPlugin_GHS3D_i::~GHS3DPlugin_GHS3D_i
+ *
+ *  Destructor
+ */
+//=============================================================================
+
+bool GHS3DPlugin_GHS3D_i::importGMFMesh(const char* theGMFFileName)
+{
+  MESSAGE( "GHS3DPlugin_GHS3D_i::importGMFMesh" );
+
+  SMESH::SMESH_Mesh_ptr theMesh = SMESH_Gen_i::GetSMESHGen()->CreateEmptyMesh();
+  SMESH_Gen_i::GetSMESHGen()->RemoveLastFromPythonScript(SMESH_Gen_i::GetSMESHGen()->GetCurrentStudy()->StudyId());
+  SALOMEDS::SObject_ptr theSMesh = SMESH_Gen_i::GetSMESHGen()->ObjectToSObject(SMESH_Gen_i::GetSMESHGen()->GetCurrentStudy(), theMesh);
+#ifdef WINNT
+#define SEP '\\'
+#else
+#define SEP '/'
+#endif
+  string strFileName (theGMFFileName);
+  strFileName = strFileName.substr(strFileName.rfind(SEP)+1);
+  strFileName.erase(strFileName.rfind('.'));
+  SMESH_Gen_i::GetSMESHGen()->SetName(theSMesh, strFileName.c_str());
+  SMESH_Mesh_i* meshServant = dynamic_cast<SMESH_Mesh_i*>( SMESH_Gen_i::GetSMESHGen()->GetServant( theMesh ).in() );
+  ASSERT( meshServant );
+  if ( meshServant ) {
+    bool res = GetImpl()->importGMFMesh(theGMFFileName, meshServant->GetImpl());
+    SMESH::TPythonDump() << "isDone = " << _this() << ".importGMFMesh( \"" << theGMFFileName << "\")";
+    return res;
+  }
+  return false;
+}
index e61fb45cdd7eef386e38189af9b23a355c2b9d59..9de897f9e2dffce7adfca13b4ae577c248568799 100644 (file)
@@ -49,6 +49,8 @@ public:
  
   // Get implementation
   ::GHS3DPlugin_GHS3D* GetImpl();
+
+  virtual bool importGMFMesh(const char* theGMFFileName);
 };
 
 #endif
index bdfa2b9fe502ad6fbd6b9436eabcc0c49d513f83..7498235c9639956d444167dae0d852bc79b90445 100644 (file)
@@ -24,6 +24,8 @@
 //=============================================================================
 //
 #include "GHS3DPlugin_Hypothesis.hxx"
+#include <SMESH_ProxyMesh.hxx>
+#include <StdMeshers_QuadToTriaAdaptor.hxx>
 
 #include <TCollection_AsciiString.hxx>
 
@@ -49,10 +51,19 @@ GHS3DPlugin_Hypothesis::GHS3DPlugin_Hypothesis(int hypId, int studyId, SMESH_Gen
   myToUseBoundaryRecoveryVersion(DefaultToUseBoundaryRecoveryVersion()),
   myToUseFemCorrection(DefaultToUseFEMCorrection()),
   myToRemoveCentralPoint(DefaultToRemoveCentralPoint()),
-  myEnforcedVertices(DefaultEnforcedVertices())
+  myEnforcedVertices(DefaultEnforcedVertices()),
+  _enfNodes(DefaultIDSortedNodeSet()),
+  _enfEdges(DefaultIDSortedElemSet()),
+  _enfTriangles(DefaultIDSortedElemSet()),
+  _enfQuadrangles(DefaultIDSortedElemSet())
 {
   _name = "GHS3D_Parameters";
   _param_algo_dim = 3;
+  _edgeID2nodeIDMap.clear();
+  _triID2nodeIDMap.clear();
+  _quadID2nodeIDMap.clear();
+  _nodeIDToSizeMap.clear();
+  _elementIDToSizeMap.clear();
 }
 
 //=======================================================================
@@ -328,6 +339,181 @@ void GHS3DPlugin_Hypothesis::SetEnforcedVertex(double x, double y, double z, dou
   NotifySubMeshesHypothesisModification();
 }
 
+
+//=======================================================================
+//function : SetEnforcedMesh
+//=======================================================================
+void GHS3DPlugin_Hypothesis::SetEnforcedMesh(SMESH_Mesh& theMesh, SMESH::ElementType elementType, double size)
+{
+  TIDSortedElemSet theElemSet;
+  SMDS_ElemIteratorPtr eIt;
+/*
+  if ((elementType == SMESH::FACE) && (theMesh.NbQuadrangles() > 0)) {
+    SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
+
+    StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor;
+    aQuad2Trias->Compute( theMesh );
+    proxyMesh.reset(aQuad2Trias );
+
+//    std::cout << "proxyMesh->NbFaces(): " << proxyMesh->NbFaces() << std::endl;
+//    eIt = proxyMesh->GetFaces();
+//    if (eIt)
+//      while ( eIt->more() )
+//        theElemSet.insert( eIt->next() );
+//    else {
+//    std::cout << "********************** eIt == 0 *****************" << std::endl;
+    eIt = theMesh.GetMeshDS()->elementsIterator(SMDSAbs_ElementType(elementType));
+    while ( eIt->more() ) {
+      const SMDS_MeshElement* elem = eIt->next();
+        theElemSet.insert( elem );
+    }
+  }
+
+  else
+  {
+  */
+    eIt = theMesh.GetMeshDS()->elementsIterator(SMDSAbs_ElementType(elementType));
+    while ( eIt->more() )
+      theElemSet.insert( eIt->next() );
+/*
+  }
+*/
+  MESSAGE("Add "<<theElemSet.size()<<" types["<<elementType<<"] from source mesh");
+
+  SetEnforcedElements( theElemSet, elementType, size);
+
+}
+
+//=======================================================================
+//function : SetEnforcedElements
+//=======================================================================
+void GHS3DPlugin_Hypothesis::SetEnforcedElements(TIDSortedElemSet theElemSet, SMESH::ElementType elementType, double size)
+{
+  MESSAGE("GHS3DPlugin_Hypothesis::SetEnforcedElements");
+  TIDSortedElemSet::const_iterator it = theElemSet.begin();
+  const SMDS_MeshElement* elem;
+  for (;it != theElemSet.end();++it)
+  {
+    elem = (*it);
+//     MESSAGE("Element ID: " << (*it)->GetID());
+    const SMDS_MeshNode* node = dynamic_cast<const SMDS_MeshNode*>(elem);
+    switch (elementType) {
+      case SMESH::NODE:
+        if (node) {
+//           MESSAGE("This is a node");
+          _enfNodes.insert(node);
+          _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+        }
+        else {
+//           MESSAGE("This is an element");
+          _enfNodes.insert(elem->begin_nodes(),elem->end_nodes());
+          SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+          for (;nodeIt->more();)
+            _nodeIDToSizeMap.insert(make_pair(nodeIt->next()->GetID(), size));
+        }
+        NotifySubMeshesHypothesisModification();
+        break;
+      case SMESH::EDGE:
+        if (node) {
+//           MESSAGE("This is a node");
+        }
+        else {
+//           MESSAGE("This is an element");
+          if (elem->GetType() == SMDSAbs_Edge) {
+            _enfEdges.insert(elem);
+//             _enfNodes.insert(elem->begin_nodes(),elem->end_nodes());
+            _elementIDToSizeMap.insert(make_pair(elem->GetID(), size));
+            SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+            for (;nodeIt->more();) {
+              node = dynamic_cast<const SMDS_MeshNode*>(nodeIt->next());
+              _edgeID2nodeIDMap[elem->GetID()].push_back(node->GetID());
+              _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+            }
+          }
+          else if (elem->GetType() > SMDSAbs_Edge) {
+            SMDS_ElemIteratorPtr it = elem->edgesIterator();
+            for (;it->more();) {
+              const SMDS_MeshElement* anEdge = it->next();
+              _enfEdges.insert(anEdge);
+//               _enfNodes.insert(anEdge->begin_nodes(),anEdge->end_nodes());
+              _elementIDToSizeMap.insert(make_pair(anEdge->GetID(), size));
+              SMDS_ElemIteratorPtr nodeIt = anEdge->nodesIterator();
+              for (;nodeIt->more();) {
+                node = dynamic_cast<const SMDS_MeshNode*>(nodeIt->next());
+                _edgeID2nodeIDMap[anEdge->GetID()].push_back(node->GetID());
+                _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+              }
+            }
+          }
+        }
+        NotifySubMeshesHypothesisModification();
+        break;
+      case SMESH::FACE:
+        if (node) {
+//           MESSAGE("This is a node");
+        }
+        else {
+//           MESSAGE("This is an element");
+          if (elem->GetType() == SMDSAbs_Face)
+          {
+            if (elem->NbNodes() == 3) {
+              _enfTriangles.insert(elem);
+//               _enfNodes.insert(elem->begin_nodes(),elem->end_nodes());
+              _elementIDToSizeMap.insert(make_pair(elem->GetID(), size));
+              SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+              for (;nodeIt->more();) {
+                node = dynamic_cast<const SMDS_MeshNode*>(nodeIt->next());
+                _triID2nodeIDMap[elem->GetID()].push_back(node->GetID());
+                _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+              }
+            }
+            else if (elem->NbNodes() == 4) {
+              _enfQuadrangles.insert(elem);
+//               _enfNodes.insert(elem->begin_nodes(),elem->end_nodes());
+              _elementIDToSizeMap.insert(make_pair(elem->GetID(), size));
+              SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+              for (;nodeIt->more();) {
+                node = dynamic_cast<const SMDS_MeshNode*>(nodeIt->next());
+                _quadID2nodeIDMap[elem->GetID()].push_back(node->GetID());
+                _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+              }
+            }
+          }
+          else if (elem->GetType() > SMDSAbs_Face) {
+            SMDS_ElemIteratorPtr it = elem->facesIterator();
+            for (;it->more();) {
+              const SMDS_MeshElement* aFace = it->next();
+              if (aFace->NbNodes() == 3) {
+                _enfTriangles.insert(aFace);
+//                 _enfNodes.insert(aFace->begin_nodes(),aFace->end_nodes());
+                _elementIDToSizeMap.insert(make_pair(aFace->GetID(), size));
+                SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+                for (;nodeIt->more();) {
+                  node = dynamic_cast<const SMDS_MeshNode*>(nodeIt->next());
+                  _triID2nodeIDMap[aFace->GetID()].push_back(node->GetID());
+                  _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+                }
+              }
+              else if (aFace->NbNodes() == 4) {
+                _enfQuadrangles.insert(aFace);
+//                 _enfNodes.insert(aFace->begin_nodes(),aFace->end_nodes());
+                _elementIDToSizeMap.insert(make_pair(aFace->GetID(), size));
+                SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+                for (;nodeIt->more();) {
+                  node = dynamic_cast<const SMDS_MeshNode*>(nodeIt->next());
+                  _quadID2nodeIDMap[aFace->GetID()].push_back(node->GetID());
+                  _nodeIDToSizeMap.insert(make_pair(node->GetID(), size));
+                }
+              }
+            }
+          }
+        }
+        NotifySubMeshesHypothesisModification();
+        break;
+    };
+  }
+}
+
 //=======================================================================
 //function : GetEnforcedVertex
 //=======================================================================
@@ -377,6 +563,24 @@ void GHS3DPlugin_Hypothesis::ClearEnforcedVertices()
     NotifySubMeshesHypothesisModification();
 }
 
+//=======================================================================
+//function : ClearEnforcedMeshes
+//=======================================================================
+void GHS3DPlugin_Hypothesis::ClearEnforcedMeshes()
+{
+   _enfNodes.clear();
+   _enfEdges.clear();
+   _enfTriangles.clear();
+   _enfQuadrangles.clear();
+//    _edgeID2nodeIDMap.clear();
+//    _triID2nodeIDMap.clear();
+//    _quadID2nodeIDMap.clear();
+//    _nodeIDToSizeMap.clear();
+//    _elementIDToSizeMap.clear();
+   NotifySubMeshesHypothesisModification();
+}
+
+
 //=======================================================================
 //function : DefaultMeshHoles
 //=======================================================================
@@ -524,6 +728,15 @@ GHS3DPlugin_Hypothesis::TEnforcedVertexValues GHS3DPlugin_Hypothesis::DefaultEnf
   return GHS3DPlugin_Hypothesis::TEnforcedVertexValues();
 }
 
+TIDSortedNodeSet GHS3DPlugin_Hypothesis::DefaultIDSortedNodeSet()
+{
+  return TIDSortedNodeSet();
+}
+
+TIDSortedElemSet GHS3DPlugin_Hypothesis::DefaultIDSortedElemSet()
+{
+  return TIDSortedElemSet();
+}
 
 //=======================================================================
 //function : SaveTo
@@ -744,7 +957,7 @@ bool GHS3DPlugin_Hypothesis::SetParametersByDefaults(const TDefaults&  /*dflts*/
 std::string GHS3DPlugin_Hypothesis::CommandToRun(const GHS3DPlugin_Hypothesis* hyp,
                                             const bool                    hasShapeToMesh)
 {
-  TCollection_AsciiString cmd( "ghs3d" );
+  TCollection_AsciiString cmd( "ghs3d" ); // ghs3dV4.1_32bits (to permit the .mesh output)
   // check if any option is overridden by hyp->myTextOption
   bool m   = hyp ? ( hyp->myTextOption.find("-m")  == std::string::npos ) : true;
   bool M   = hyp ? ( hyp->myTextOption.find("-M")  == std::string::npos ) : true;
@@ -881,3 +1094,47 @@ GHS3DPlugin_Hypothesis::TEnforcedVertexValues GHS3DPlugin_Hypothesis::GetEnforce
     return hyp ? hyp->_GetEnforcedVertices():DefaultEnforcedVertices();
 }
 
+TIDSortedNodeSet GHS3DPlugin_Hypothesis::GetEnforcedNodes(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetEnforcedNodes():DefaultIDSortedNodeSet();
+}
+
+TIDSortedElemSet GHS3DPlugin_Hypothesis::GetEnforcedEdges(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetEnforcedEdges():DefaultIDSortedElemSet();
+}
+
+TIDSortedElemSet GHS3DPlugin_Hypothesis::GetEnforcedTriangles(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetEnforcedTriangles():DefaultIDSortedElemSet();
+}
+
+TIDSortedElemSet GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetEnforcedQuadrangles():DefaultIDSortedElemSet();
+}
+
+GHS3DPlugin_Hypothesis::TElemID2NodeIDMap GHS3DPlugin_Hypothesis::GetEdgeID2NodeIDMap(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetEdgeID2NodeIDMap(): GHS3DPlugin_Hypothesis::TElemID2NodeIDMap();
+}
+
+GHS3DPlugin_Hypothesis::TElemID2NodeIDMap GHS3DPlugin_Hypothesis::GetTri2NodeMap(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetTri2NodeMap(): GHS3DPlugin_Hypothesis::TElemID2NodeIDMap();
+}
+
+GHS3DPlugin_Hypothesis::TElemID2NodeIDMap GHS3DPlugin_Hypothesis::GetQuad2NodeMap(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetQuad2NodeMap(): GHS3DPlugin_Hypothesis::TElemID2NodeIDMap();
+}
+
+GHS3DPlugin_Hypothesis::TID2SizeMap GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetNodeIDToSizeMap(): GHS3DPlugin_Hypothesis::TID2SizeMap();
+}
+
+GHS3DPlugin_Hypothesis::TID2SizeMap GHS3DPlugin_Hypothesis::GetElementIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp)
+{
+    return hyp ? hyp->_GetElementIDToSizeMap(): GHS3DPlugin_Hypothesis::TID2SizeMap();
+}
index 28fe2bb64b201c5d9c2c4e9a175bb3870b152353..ac546762dc406d49a560e7553529d7eee41ed733 100644 (file)
 
 #include "GHS3DPlugin_Defs.hxx"
 
-#include <SMESH_Hypothesis.hxx>
-#include <utilities.h>
+#include "SMESH_Hypothesis.hxx"
+#include "SMESH_Mesh_i.hxx"
+#include "SMESH_Gen_i.hxx"
+#include "SMESH_TypeDefs.hxx"
+#include "utilities.h"
 
 #include <stdexcept>
 #include <map>
@@ -117,6 +120,39 @@ public:
   /*!
    * To set an enforced vertex
    */
+//   struct TEnforcedNode {
+//     std::vector<double> coords;
+//     double size;
+//     std::string geomEntry;
+//     std::string groupName;
+//   };
+//   
+//   struct CompareEnfNodes {
+//     bool operator () (const TEnforcedNode* e1, const TEnforcedNode* e2) const {
+//       if (e1 && e2) {
+//         if (e1->coords.size() && e2->coords.size())
+//           return (e1->coords < e2->coords);
+//         else
+//           return (e1->geomEntry < e2->geomEntry);
+//       }
+//       return false;
+//     }
+//   };
+//   typedef std::set< TEnforcedNode*, CompareEnfNodes > TEnforcedNodeList;
+//   // Map Coords / Enforced node
+//   typedef std::map< std::vector<double>, TEnforcedNode* > TCoordsEnfNodeMap;
+//   // Map geom entry / Enforced ndoe
+//   typedef std::map< std::string, TEnforcedNode* > TGeomEntryEnfNodeMap;
+//   
+//   
+//   struct TEnforcedEdge {
+//     long ID;
+//     long node1;
+//     long node2;
+//     std::string groupName;
+//   };
+  
+  
   typedef std::map<std::vector<double>,double> TEnforcedVertexValues;
   void SetEnforcedVertex(double x, double y, double z, double size);
   double GetEnforcedVertex(double x, double y, double z) throw (std::invalid_argument);
@@ -124,6 +160,24 @@ public:
   const TEnforcedVertexValues _GetEnforcedVertices() const { return myEnforcedVertices; }
   void ClearEnforcedVertices();
 
+  /*!
+   * To set enforced elements
+   */
+  void SetEnforcedMesh(SMESH_Mesh& theMesh, SMESH::ElementType elementType, double size);
+  void SetEnforcedElements(TIDSortedElemSet theElemSet, SMESH::ElementType elementType, double size);
+  void ClearEnforcedMeshes();
+  const TIDSortedNodeSet _GetEnforcedNodes() const { return _enfNodes; }
+  const TIDSortedElemSet _GetEnforcedEdges() const { return _enfEdges; }
+  const TIDSortedElemSet _GetEnforcedTriangles() const { return _enfTriangles; }
+  const TIDSortedElemSet _GetEnforcedQuadrangles() const { return _enfQuadrangles; }
+  typedef std::map<int,std::vector<int> > TElemID2NodeIDMap;
+  const TElemID2NodeIDMap _GetEdgeID2NodeIDMap() const {return _edgeID2nodeIDMap; }
+  const TElemID2NodeIDMap _GetTri2NodeMap() const {return _triID2nodeIDMap; }
+  const TElemID2NodeIDMap _GetQuad2NodeMap() const {return _quadID2nodeIDMap; }
+  typedef std::map<int,double> TID2SizeMap;
+  const TID2SizeMap _GetNodeIDToSizeMap() const {return _nodeIDToSizeMap; }
+  const TID2SizeMap _GetElementIDToSizeMap() const {return _elementIDToSizeMap; }
+
   static bool   DefaultMeshHoles();
   static short  DefaultMaximumMemory();
   static short  DefaultInitialMemory();
@@ -136,6 +190,8 @@ public:
   static bool   DefaultToUseFEMCorrection();
   static bool   DefaultToRemoveCentralPoint();
   static TEnforcedVertexValues DefaultEnforcedVertices();
+  static TIDSortedNodeSet DefaultIDSortedNodeSet();
+  static TIDSortedElemSet DefaultIDSortedElemSet();
 
   /*!
    * \brief Return command to run ghs3d mesher excluding file prefix (-f)
@@ -150,7 +206,16 @@ public:
    * \brief Return the enforced vertices
    */
   static TEnforcedVertexValues GetEnforcedVertices(const GHS3DPlugin_Hypothesis* hyp);
-
+  static TIDSortedNodeSet GetEnforcedNodes(const GHS3DPlugin_Hypothesis* hyp);
+  static TIDSortedElemSet GetEnforcedEdges(const GHS3DPlugin_Hypothesis* hyp);
+  static TIDSortedElemSet GetEnforcedTriangles(const GHS3DPlugin_Hypothesis* hyp);
+  static TIDSortedElemSet GetEnforcedQuadrangles(const GHS3DPlugin_Hypothesis* hyp);
+  static TElemID2NodeIDMap GetEdgeID2NodeIDMap(const GHS3DPlugin_Hypothesis* hyp);
+  static TElemID2NodeIDMap GetTri2NodeMap(const GHS3DPlugin_Hypothesis* hyp);
+  static TElemID2NodeIDMap GetQuad2NodeMap(const GHS3DPlugin_Hypothesis* hyp);
+  static TID2SizeMap GetNodeIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp);
+  static TID2SizeMap GetElementIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp);
+  
   // Persistence
   virtual std::ostream & SaveTo(std::ostream & save);
   virtual std::istream & LoadFrom(std::istream & load);
@@ -182,7 +247,15 @@ private:
   bool   myToRemoveCentralPoint;
   std::string myTextOption;
   TEnforcedVertexValues myEnforcedVertices;
-  
+  TIDSortedNodeSet _enfNodes;
+  TIDSortedElemSet _enfEdges;
+  TIDSortedElemSet _enfTriangles;
+  TIDSortedElemSet _enfQuadrangles;
+  TElemID2NodeIDMap _edgeID2nodeIDMap;
+  TElemID2NodeIDMap _triID2nodeIDMap;
+  TElemID2NodeIDMap _quadID2nodeIDMap;
+  TID2SizeMap _nodeIDToSizeMap;
+  TID2SizeMap _elementIDToSizeMap;
 };
 
 
index ff3953e3d1fda035db4dc22733064dde56c62a16..63952ba28f41805f3a6ed2aa454802407de0d5bd 100644 (file)
 //
 #include "GHS3DPlugin_Hypothesis_i.hxx"
 
-#include <SMESH_Gen.hxx>
-#include <SMESH_PythonDump.hxx>
-
-#include <Utils_CorbaException.hxx>
-#include <utilities.h>
-#include <SMESH_Mesh_i.hxx>
-
+#include "SMESH_Gen.hxx"
+#include "SMESH_PythonDump.hxx"
+//#include "SMESH_Mesh.hxx"
+//#include "SMESH_ProxyMesh.hxx"
+//#include <StdMeshers_QuadToTriaAdaptor.hxx>
+
+#include "Utils_CorbaException.hxx"
+#include "utilities.h"
+#include "SMESH_Mesh_i.hxx"
+#include "SMESH_Group_i.hxx"
+#include "SMESH_Gen_i.hxx"
+#include "SMESH_TypeDefs.hxx"
+
+#ifndef GHS3D_VERSION
+#define GHS3D_VERSION 41
+#endif
 //=======================================================================
 //function : GHS3DPlugin_Hypothesis_i
 //=======================================================================
@@ -418,7 +427,7 @@ void GHS3DPlugin_Hypothesis_i::RemoveEnforcedVertex(CORBA::Double x, CORBA::Doub
     ExDescription.text = ex.what();
     ExDescription.type = SALOME::BAD_PARAM;
     ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::RemoveEnforcedVertex(x,y,z)";
-    ExDescription.lineNumber = 0;
+    ExDescription.lineNumber = 408;
     throw SALOME::SALOME_Exception(ExDescription);
   }
   catch (SALOME_Exception& ex) {
@@ -434,8 +443,195 @@ void GHS3DPlugin_Hypothesis_i::ClearEnforcedVertices()
 {
   ASSERT(myBaseImpl);
   this->GetImpl()->ClearEnforcedVertices();
+  SMESH::TPythonDump () << _this() << ".ClearEnforcedVertices() ";
+}
+
+//=======================================================================
+//function : ClearEnforcedMeshes
+//=======================================================================
+
+void GHS3DPlugin_Hypothesis_i::ClearEnforcedMeshes()
+{
+  ASSERT(myBaseImpl);
+  this->GetImpl()->ClearEnforcedMeshes();
+  SMESH::TPythonDump () << _this() << ".ClearEnforcedMeshes() ";
+}
+
+/*!
+ * \brief Adds enforced elements of type elementType using another mesh/sub-mesh/mesh group theSource.
+ */
+void GHS3DPlugin_Hypothesis_i::SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType theType)
+  throw (SALOME::SALOME_Exception)
+{
+#if GHS3D_VERSION >= 42
+  _SetEnforcedMesh(theSource, theType, -1.0);
+  SMESH_Mesh_i* theMesh_i = SMESH::DownCast<SMESH_Mesh_i*>( theSource);
+  SMESH_Group_i* theGroup_i = SMESH::DownCast<SMESH_Group_i*>( theSource);
+  if (theGroup_i)
+  {
+    SMESH::TPythonDump () << _this() << ".SetEnforcedMesh( " 
+                          << theSource << ", " << theType << " )";
+  }
+  else if (theMesh_i)
+  {
+    SMESH::TPythonDump () << _this() << ".SetEnforcedMesh( " 
+                          << theSource << ".GetMesh(), " << theType << " )";
+  }
+#else
+  SALOME::ExceptionStruct ExDescription;
+  ExDescription.text = "Bad version of GHS3D. It must >= 4.2.";
+  ExDescription.type = SALOME::BAD_PARAM;
+  ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::SetEnforcedMesh(theSource, theType)";
+  ExDescription.lineNumber = 463;
+  throw SALOME::SALOME_Exception(ExDescription);
+#endif
+}
+
+/*!
+ * \brief Adds enforced elements of type elementType using another mesh/sub-mesh/mesh group theSource and a size.
+ */
+void GHS3DPlugin_Hypothesis_i::SetEnforcedMeshSize(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType theType, double theSize)
+  throw (SALOME::SALOME_Exception)
+{
+  if (theSize <= 0) {
+    SALOME::ExceptionStruct ExDescription;
+    ExDescription.text = "Size cannot be negative";
+    ExDescription.type = SALOME::BAD_PARAM;
+    ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::SetEnforcedMeshSize(theSource, theType)";
+    ExDescription.lineNumber = 475;
+    throw SALOME::SALOME_Exception(ExDescription);
+  }
+  
+  _SetEnforcedMesh(theSource, theType, theSize);
+  SMESH_Mesh_i* theMesh_i = SMESH::DownCast<SMESH_Mesh_i*>( theSource);
+  SMESH_Group_i* theGroup_i = SMESH::DownCast<SMESH_Group_i*>( theSource);
+  if (theGroup_i)
+  {
+    SMESH::TPythonDump () << _this() << ".SetEnforcedMeshSize( " 
+                          << theSource << ", " << theType << ", " << theSize << " )";
+  }
+  else if (theMesh_i)
+  {
+    SMESH::TPythonDump () << _this() << ".SetEnforcedMeshSize( " 
+                          << theSource << ".GetMesh(), " << theType << ", " << theSize << " )";
+  }
 }
 
+void GHS3DPlugin_Hypothesis_i::_SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType theType, double theSize)
+  throw (SALOME::SALOME_Exception)
+{
+  ASSERT(myBaseImpl);
+  
+  if (CORBA::is_nil( theSource ))
+  {
+    SALOME::ExceptionStruct ExDescription;
+    ExDescription.text = "The source mesh CORBA object is NULL";
+    ExDescription.type = SALOME::BAD_PARAM;
+    ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)";
+    ExDescription.lineNumber = 502;
+    throw SALOME::SALOME_Exception(ExDescription);
+  }
+  
+  if ((theType != SMESH::NODE) && (theType != SMESH::EDGE) && (theType != SMESH::FACE))
+  {
+    SALOME::ExceptionStruct ExDescription;
+    ExDescription.text = "Bad elementType";
+    ExDescription.type = SALOME::BAD_PARAM;
+    ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)";
+    ExDescription.lineNumber = 502;
+    throw SALOME::SALOME_Exception(ExDescription);
+  }
+  
+  SMESH::array_of_ElementType_var types = theSource->GetTypes();
+//   MESSAGE("Required type is "<<theType);
+//   MESSAGE("Available types:");
+//   for (int i=0;i<types->length();i++){MESSAGE(types[i]);}
+  if ( types->length() >= 1 && types[types->length()-1] <  theType)
+  {
+    SALOME::ExceptionStruct ExDescription;
+    ExDescription.text = "The source mesh has bad type";
+    ExDescription.type = SALOME::BAD_PARAM;
+    ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)";
+    ExDescription.lineNumber = 502;
+    throw SALOME::SALOME_Exception(ExDescription);
+  }
+  
+  SMESHDS_Mesh* theMeshDS;
+  SMESH_Mesh_i* anImplPtr = SMESH::DownCast<SMESH_Mesh_i*>(theSource->GetMesh());
+  if (anImplPtr)
+    theMeshDS = anImplPtr->GetImpl().GetMeshDS();
+  else
+    return;
+  
+  SMESH_Mesh_i* theMesh_i = SMESH::DownCast<SMESH_Mesh_i*>( theSource);
+  SMESH_Group_i* theGroup_i = SMESH::DownCast<SMESH_Group_i*>( theSource);
+  TIDSortedElemSet theElemSet;
+
+  if (theMesh_i)
+  {
+    try {
+      this->GetImpl()->SetEnforcedMesh(anImplPtr->GetImpl(), theType, theSize);
+    }
+    catch (const std::invalid_argument& ex) {
+      SALOME::ExceptionStruct ExDescription;
+      ExDescription.text = ex.what();
+      ExDescription.type = SALOME::BAD_PARAM;
+      ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)";
+      ExDescription.lineNumber = 502;
+      throw SALOME::SALOME_Exception(ExDescription);
+    }
+    catch (SALOME_Exception& ex) {
+      THROW_SALOME_CORBA_EXCEPTION( ex.what() ,SALOME::BAD_PARAM );
+    }
+//
+//      SMESH::long_array_var anIDs = theMesh_i->GetElementsByType(theType);
+//      if ( anIDs->length() == 0 ){MESSAGE("The source mesh is empty");}
+//      for (int i=0; i<anIDs->length(); i++) {
+//        CORBA::Long ind = anIDs[i];
+//        const SMDS_MeshElement * elem = theMeshDS->FindElement(ind);
+//        if (elem)
+//          theElemSet.insert( elem );
+//      }
+////    }
+//    MESSAGE("Add "<<theElemSet.size()<<" types["<<theType<<"] from source mesh");
+  }
+  else if (theGroup_i && types->length() == 1 && types[0] == theType)
+  {
+    SMESH::long_array_var anIDs = theGroup_i->GetListOfID();
+    if ( anIDs->length() == 0 ){MESSAGE("The source group is empty");}
+    for (int i=0; i<anIDs->length(); i++) {
+      CORBA::Long ind = anIDs[i];
+      if (theType == SMESH::NODE)
+      {
+        const SMDS_MeshNode * node = theMeshDS->FindNode(ind);
+        if (node)
+          theElemSet.insert( node );
+      }
+      else
+      {
+        const SMDS_MeshElement * elem = theMeshDS->FindElement(ind);
+        if (elem)
+          theElemSet.insert( elem );
+      }
+    }
+    MESSAGE("Add "<<theElemSet.size()<<" types["<<theType<<"] from source group "<< theGroup_i->GetName());
+
+    try {
+      this->GetImpl()->SetEnforcedElements(theElemSet, theType, theSize);
+    }
+    catch (const std::invalid_argument& ex) {
+      SALOME::ExceptionStruct ExDescription;
+      ExDescription.text = ex.what();
+      ExDescription.type = SALOME::BAD_PARAM;
+      ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)";
+      ExDescription.lineNumber = 502;
+      throw SALOME::SALOME_Exception(ExDescription);
+    }
+    catch (SALOME_Exception& ex) {
+      THROW_SALOME_CORBA_EXCEPTION( ex.what() ,SALOME::BAD_PARAM );
+    }
+  }
+}
 //=============================================================================
 /*!
  *  Get implementation
index 0cf9a26606649d196b30b49cc5c604008509775f..d9b59414e8b403b282902b4a87cefac188bcfda8 100644 (file)
@@ -31,6 +31,7 @@
 #include CORBA_SERVER_HEADER(GHS3DPlugin_Algorithm)
 
 #include "SMESH_Hypothesis_i.hxx"
+#include "SMESH_Mesh_i.hxx"
 #include "GHS3DPlugin_Hypothesis.hxx"
 
 class SMESH_Gen;
@@ -129,6 +130,13 @@ class GHS3DPLUGIN_EXPORT GHS3DPlugin_Hypothesis_i:
   void RemoveEnforcedVertex(CORBA::Double x, CORBA::Double y, CORBA::Double z) throw (SALOME::SALOME_Exception);
   GHS3DPlugin::GHS3DEnforcedVertexList* GetEnforcedVertices();
   void ClearEnforcedVertices();
+  /*!
+   * To set an enforced mesh
+   */
+  void SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType elementType) throw (SALOME::SALOME_Exception);
+  void SetEnforcedMeshSize(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType elementType, double size) throw (SALOME::SALOME_Exception);
+  void _SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType elementType, double size) throw (SALOME::SALOME_Exception);
+  void ClearEnforcedMeshes();
 
   // Get implementation
   ::GHS3DPlugin_Hypothesis* GetImpl();
index 28751997741112da0836d41051ec9bceccee53e8..4e7dda2be4569a59a72da0550cfec0723859a897 100644 (file)
@@ -27,6 +27,7 @@ include $(top_srcdir)/adm_local/unix/make_common_starter.am
 
 # header files 
 salomeinclude_HEADERS = \
+       libmesh5.h \
        GHS3DPlugin_Defs.hxx \
        GHS3DPlugin_GHS3D.hxx \
        GHS3DPlugin_GHS3D_i.hxx \
@@ -37,6 +38,7 @@ salomeinclude_HEADERS = \
 lib_LTLIBRARIES = libGHS3DEngine.la
 
 dist_libGHS3DEngine_la_SOURCES = \
+  libmesh5.c \
   GHS3DPlugin_GHS3D.cxx \
   GHS3DPlugin_GHS3D_i.cxx \
   GHS3DPlugin_i.cxx \
diff --git a/src/GHS3DPlugin/libmesh5.c b/src/GHS3DPlugin/libmesh5.c
new file mode 100644 (file)
index 0000000..4277b97
--- /dev/null
@@ -0,0 +1,1192 @@
+
+
+/*----------------------------------------------------------*/
+/*                                                                                                                     */
+/*                                             LIBMESH V 5.45                                          */
+/*                                                                                                                     */
+/*----------------------------------------------------------*/
+/*                                                                                                                     */
+/*     Description:            handle .meshb file format I/O           */
+/*     Author:                         Loic MARECHAL                                           */
+/*     Creation date:          feb 16 2007                                                     */
+/*     Last modification:      feb 08 2011                                                     */
+/*                                                                                                                     */
+/*----------------------------------------------------------*/
+
+
+/*----------------------------------------------------------*/
+/* Includes                                                                                                    */
+/*----------------------------------------------------------*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include <ctype.h>
+#include "libmesh5.h"
+
+
+/*----------------------------------------------------------*/
+/* Defines                                                                                                     */
+/*----------------------------------------------------------*/
+
+#define Asc 1
+#define Bin 2
+#define MshFil 4
+#define SolFil 8
+#define MaxMsh 100
+#define InfKwd 1
+#define RegKwd 2
+#define SolKwd 3
+#define WrdSiz 4
+#define BufSiz 10000
+
+
+/*----------------------------------------------------------*/
+/* Structures                                                                                          */
+/*----------------------------------------------------------*/
+
+typedef struct
+{
+       int typ, SolSiz, NmbWrd, NmbLin, NmbTyp, TypTab[ GmfMaxTyp ];
+       long pos;
+       char fmt[ GmfMaxTyp*9 ];
+}KwdSct;
+
+typedef struct
+{
+       int dim, ver, mod, typ, cod, pos;
+       long NexKwdPos;
+       KwdSct KwdTab[ GmfMaxKwd + 1 ];
+       FILE *hdl;
+       int *IntBuf;
+       float *FltBuf;
+       unsigned char *buf;
+       char FilNam[ GmfStrSiz ];
+       double DblBuf[1000/8];
+       unsigned char blk[ BufSiz + 1000 ];
+}GmfMshSct;
+
+
+/*----------------------------------------------------------*/
+/* Global variables                                                                                    */
+/*----------------------------------------------------------*/
+
+int GmfIniFlg=0;
+GmfMshSct *GmfMshTab[ MaxMsh + 1 ];
+char *GmfKwdFmt[ GmfMaxKwd + 1 ][4] = 
+{      {"Reserved", "", "", ""},
+       {"MeshVersionFormatted", "", "", "i"},
+       {"Reserved", "", "", ""},
+       {"Dimension", "", "", "i"},
+       {"Vertices", "Vertex", "i", "dri"},
+       {"Edges", "Edge", "i", "iii"},
+       {"Triangles", "Triangle", "i", "iiii"},
+       {"Quadrilaterals", "Quadrilateral", "i", "iiiii"},
+       {"Tetrahedra", "Tetrahedron", "i", "iiiii"},
+       {"Prisms", "Prism", "i", "iiiiiii"},
+       {"Hexahedra", "Hexahedron", "i", "iiiiiiiii"},
+       {"IterationsAll", "IterationAll","","i"},
+       {"TimesAll", "TimeAll","","r"},                                 
+       {"Corners", "Corner", "i", "i"},
+       {"Ridges", "Ridge", "i", "i"},
+       {"RequiredVertices", "RequiredVertex", "i", "i"},
+       {"RequiredEdges", "RequiredEdge", "i", "i"},
+       {"RequiredTriangles", "RequiredTriangle", "i", "i"},
+       {"RequiredQuadrilaterals", "RequiredQuadrilateral", "i", "i"},
+       {"TangentAtEdgeVertices", "TangentAtEdgeVertex", "i", "iii"},
+       {"NormalAtVertices", "NormalAtVertex", "i", "ii"},
+       {"NormalAtTriangleVertices", "NormalAtTriangleVertex", "i", "iii"},
+       {"NormalAtQuadrilateralVertices", "NormalAtQuadrilateralVertex", "i", "iiii"},
+       {"AngleOfCornerBound", "", "", "r"},
+       {"TrianglesP2", "TriangleP2", "i", "iiiiiii"},
+       {"EdgesP2", "EdgeP2", "i", "iiii"},
+       {"SolAtPyramids", "SolAtPyramid", "i", "sr"},
+       {"QuadrilateralsQ2", "QuadrilateralQ2", "i", "iiiiiiiiii"},
+       {"ISolAtPyramids", "ISolAtPyramid", "i", "iiiii"},
+       {"Reserved", "", "", ""},
+       {"TetrahedraP2", "TetrahedronP2", "i", "iiiiiiiiiii"},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"HexahedraQ2", "HexahedronQ2", "i", "iiiiiiiiiiiiiiiiiiiiiiiiiiii"},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Polyhedra", "Polyhedron", "i", "iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"},
+       {"Polygons", "Polygon", "", "iiiiiiiii"},
+       {"Reserved", "", "", ""},
+       {"Pyramids", "Pyramid", "i", "iiiiii"},
+       {"BoundingBox", "", "", "drdr"},
+       {"Body","i", "drdrdrdr"},
+       {"PrivateTable", "PrivateTable", "i", "i"},
+       {"Reserved", "", "", ""},
+       {"End", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Reserved", "", "", ""},
+       {"Tangents", "Tangent", "i", "dr"},
+       {"Normals", "Normal", "i", "dr"},
+       {"TangentAtVertices", "TangentAtVertex", "i", "ii"},
+       {"SolAtVertices", "SolAtVertex", "i", "sr"},
+       {"SolAtEdges", "SolAtEdge", "i", "sr"},
+       {"SolAtTriangles", "SolAtTriangle", "i", "sr"},
+       {"SolAtQuadrilaterals", "SolAtQuadrilateral", "i", "sr"},
+       {"SolAtTetrahedra", "SolAtTetrahedron", "i", "sr"},
+       {"SolAtPrisms", "SolAtPrism", "i", "sr"},
+       {"SolAtHexahedra", "SolAtHexahedron", "i", "sr"},
+       {"DSolAtVertices", "DSolAtVertex", "i", "sr"},
+       {"ISolAtVertices", "ISolAtVertex", "i", "i"},
+       {"ISolAtEdges", "ISolAtEdge", "i", "ii"},
+       {"ISolAtTriangles", "ISolAtTriangle", "i", "iii"},
+       {"ISolAtQuadrilaterals", "ISolAtQuadrilateral", "i", "iiii"},
+       {"ISolAtTetrahedra", "ISolAtTetrahedron", "i", "iiii"},
+       {"ISolAtPrisms", "ISolAtPrism", "i", "iiiiii"},
+       {"ISolAtHexahedra", "ISolAtHexahedron", "i", "iiiiiiii"},
+       {"Iterations", "","","i"},
+       {"Time", "","","r"},
+       {"Reserved", "","",""}
+ };
+
+
+/*----------------------------------------------------------*/
+/* Prototypes of local procedures                                                      */
+/*----------------------------------------------------------*/
+
+static void ScaWrd(GmfMshSct *, unsigned char *);
+static void ScaDblWrd(GmfMshSct *, unsigned char *);
+static void ScaBlk(GmfMshSct *, unsigned char *, int);
+static long GetPos(GmfMshSct *);
+static void RecWrd(GmfMshSct *, unsigned char *);
+static void RecDblWrd(GmfMshSct *, unsigned char *);
+static void RecBlk(GmfMshSct *, unsigned char *, int);
+static void SetPos(GmfMshSct *, long);
+static int ScaKwdTab(GmfMshSct *);
+static void ExpFmt(GmfMshSct *, int);
+static void ScaKwdHdr(GmfMshSct *, int);
+
+
+/*----------------------------------------------------------*/
+/* Open a mesh file in read or write mod                                       */
+/*----------------------------------------------------------*/
+
+int GmfOpenMesh(const char *FilNam, int mod, ...)
+{
+       int i, KwdCod, res, *PtrVer, *PtrDim, MshIdx=0;
+       char str[ GmfStrSiz ];
+       va_list VarArg;
+       GmfMshSct *msh;
+
+       if(!GmfIniFlg)
+       {
+               for(i=0;i<=MaxMsh;i++)
+                       GmfMshTab[i] = NULL;
+
+               GmfIniFlg = 1;
+       }
+
+       /*---------------------*/
+       /* MESH STRUCTURE INIT */
+       /*---------------------*/
+
+       for(i=1;i<=MaxMsh;i++)
+               if(!GmfMshTab[i])
+               {
+                       MshIdx = i;
+                       break;
+               }
+
+       if( !MshIdx || !(msh = calloc(1, sizeof(GmfMshSct))) )
+               return(0);
+
+       /* Copy the FilNam into the structure */
+
+       if(strlen(FilNam) + 7 >= GmfStrSiz)
+               return(0);
+
+       strcpy(msh->FilNam, FilNam);
+
+       /* Store the opening mod (read or write) and guess the filetype (binary or ascii) depending on the extension */
+
+       msh->mod = mod;
+       msh->buf = (void *)msh->DblBuf;
+       msh->FltBuf = (void *)msh->DblBuf;
+       msh->IntBuf = (void *)msh->DblBuf;
+
+       if(strstr(msh->FilNam, ".meshb"))
+               msh->typ |= (Bin | MshFil);
+       else if(strstr(msh->FilNam, ".mesh"))
+               msh->typ |= (Asc | MshFil);
+       else if(strstr(msh->FilNam, ".solb"))
+               msh->typ |= (Bin | SolFil);
+       else if(strstr(msh->FilNam, ".sol"))
+               msh->typ |= (Asc | SolFil);
+       else
+               return(0);
+
+       /* Open the file in the required mod and initialyse the mesh structure */
+
+       if(msh->mod == GmfRead)
+       {
+
+               /*-----------------------*/
+               /* OPEN FILE FOR READING */
+               /*-----------------------*/
+
+               va_start(VarArg, mod);
+               PtrVer = va_arg(VarArg, int *);
+               PtrDim = va_arg(VarArg, int *);
+               va_end(VarArg);
+
+               /* Create the name string and open the file */
+
+               if(!(msh->hdl = fopen(msh->FilNam, "rb")))
+                       return(0);
+
+               /* Read the endian coding tag, the mesh version and the mesh dimension (mandatory kwd) */
+
+               if(msh->typ & Bin)
+               {
+                       fread((unsigned char *)&msh->cod, WrdSiz, 1, msh->hdl);
+
+                       if( (msh->cod != 1) && (msh->cod != 16777216) )
+                               return(0);
+
+                       ScaWrd(msh, (unsigned char *)&msh->ver);
+
+                       if( (msh->ver < 1) || (msh->ver > 3) )
+                               return(0);
+
+                       if( (msh->ver == 3) && (sizeof(long) == 4) )
+                               return(0);
+
+                       ScaWrd(msh, (unsigned char *)&KwdCod);
+
+                       if(KwdCod != GmfDimension)
+                               return(0);
+
+                       GetPos(msh);
+                       ScaWrd(msh, (unsigned char *)&msh->dim);
+               }
+               else
+               {
+                       do
+                       {
+                               res = fscanf(msh->hdl, "%s", str);
+                       }while( (res != EOF) && strcmp(str, "MeshVersionFormatted") );
+
+                       if(res == EOF)
+                               return(0);
+
+                       fscanf(msh->hdl, "%d", &msh->ver);
+
+                       if( (msh->ver < 1) || (msh->ver > 3) )
+                               return(0);
+
+                       do
+                       {
+                               res = fscanf(msh->hdl, "%s", str);
+                       }while( (res != EOF) && strcmp(str, "Dimension") );
+
+                       if(res == EOF)
+                               return(0);
+
+                       fscanf(msh->hdl, "%d", &msh->dim);
+               }
+
+               if( (msh->dim != 2) && (msh->dim != 3) )
+                       return(0);
+
+               (*PtrVer) = msh->ver;
+               (*PtrDim) = msh->dim;
+
+               /*------------*/
+               /* KW READING */
+               /*------------*/
+
+               /* Read the list of kw present in the file */
+
+               if(!ScaKwdTab(msh))
+                       return(0);
+
+               GmfMshTab[ MshIdx ] = msh;
+
+               return(MshIdx);
+       }
+       else if(msh->mod == GmfWrite)
+       {
+
+               /*-----------------------*/
+               /* OPEN FILE FOR WRITING */
+               /*-----------------------*/
+
+               msh->cod = 1;
+
+               /* Check if the user provided a valid version number and dimension */
+
+               va_start(VarArg, mod);
+               msh->ver = va_arg(VarArg, int);
+               msh->dim = va_arg(VarArg, int);
+               va_end(VarArg);
+
+               if( (msh->ver < 1) || (msh->ver > 3) )
+                       return(0);
+
+               if( (msh->ver == 3) && (sizeof(long) == 4) )
+                       return(0);
+
+               if( (msh->dim != 2) && (msh->dim != 3) )
+                       return(0);
+
+               /* Create the mesh file */
+
+               if(!(msh->hdl = fopen(msh->FilNam, "wb")))
+                       return(0);
+
+               GmfMshTab[ MshIdx ] = msh;
+
+
+               /*------------*/
+               /* KW WRITING */
+               /*------------*/
+
+               /* Write the mesh version and dimension */
+
+               if(msh->typ & Asc)
+               {
+                       fprintf(msh->hdl, "%s %d\n\n", GmfKwdFmt[ GmfVersionFormatted ][0], msh->ver);
+                       fprintf(msh->hdl, "%s %d\n", GmfKwdFmt[ GmfDimension ][0], msh->dim);
+               }
+               else
+               {
+                       RecWrd(msh, (unsigned char *)&msh->cod);
+                       RecWrd(msh, (unsigned char *)&msh->ver);
+                       GmfSetKwd(MshIdx, GmfDimension, 0);
+                       RecWrd(msh, (unsigned char *)&msh->dim);
+               }
+
+               return(MshIdx);
+       }
+       else
+               return(0);
+}
+
+
+/*----------------------------------------------------------*/
+/* Close a meshfile in the right way                                           */
+/*----------------------------------------------------------*/
+
+int GmfCloseMesh(int MshIdx)
+{
+       int res = 1;
+       GmfMshSct *msh;
+
+       if( (MshIdx < 1) || (MshIdx > MaxMsh) )
+               return(0);
+
+       msh = GmfMshTab[ MshIdx ];
+       RecBlk(msh, msh->buf, 0);
+
+       /* In write down the "End" kw in write mode */
+
+  if(msh->mod == GmfWrite) {
+               if(msh->typ & Asc)
+                       fprintf(msh->hdl, "\n%s\n", GmfKwdFmt[ GmfEnd ][0]);
+               else
+                       GmfSetKwd(MshIdx, GmfEnd, 0);
+  }
+
+       /* Close the file and free the mesh structure */
+
+       if(fclose(msh->hdl))
+               res = 0;
+
+       free(msh);
+       GmfMshTab[ MshIdx ] = NULL;
+
+       return(res);
+}
+
+
+/*----------------------------------------------------------*/
+/* Read the number of lines and set the position to this kwd*/
+/*----------------------------------------------------------*/
+
+int GmfStatKwd(int MshIdx, int KwdCod, ...)
+{
+       int i, *PtrNmbTyp, *PtrSolSiz, *TypTab;
+       GmfMshSct *msh;
+       KwdSct *kwd;
+       va_list VarArg;
+
+       if( (MshIdx < 1) || (MshIdx > MaxMsh) )
+               return(0);
+
+       msh = GmfMshTab[ MshIdx ];
+
+       if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) )
+               return(0);
+
+       kwd = &msh->KwdTab[ KwdCod ];
+
+       if(!kwd->NmbLin)
+               return(0);
+
+       /* Read further arguments if this kw is a sol */
+
+       if(kwd->typ == SolKwd)
+       {
+               va_start(VarArg, KwdCod);
+
+               PtrNmbTyp = va_arg(VarArg, int *);
+               *PtrNmbTyp = kwd->NmbTyp;
+
+               PtrSolSiz = va_arg(VarArg, int *);
+               *PtrSolSiz = kwd->SolSiz;
+
+               TypTab = va_arg(VarArg, int *);
+
+               for(i=0;i<kwd->NmbTyp;i++)
+                       TypTab[i] = kwd->TypTab[i];
+
+               va_end(VarArg);
+       }
+
+       return(kwd->NmbLin);
+}
+
+
+/*----------------------------------------------------------*/
+/* Set the current file position to a given kwd                                */
+/*----------------------------------------------------------*/
+
+int GmfGotoKwd(int MshIdx, int KwdCod)
+{
+       GmfMshSct *msh;
+       KwdSct *kwd;
+
+       if( (MshIdx < 1) || (MshIdx > MaxMsh) )
+               return(0);
+
+       msh = GmfMshTab[ MshIdx ];
+
+       if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) )
+               return(0);
+
+       kwd = &msh->KwdTab[ KwdCod ];
+
+       if(!kwd->NmbLin)
+               return(0);
+
+       return(fseek(msh->hdl, kwd->pos, SEEK_SET));
+}
+
+
+/*----------------------------------------------------------*/
+/* Write the kwd and set the number of lines                           */
+/*----------------------------------------------------------*/
+
+int GmfSetKwd(int MshIdx, int KwdCod, ...)
+{
+       int i, NmbLin=0, *TypTab;
+       long CurPos;
+       va_list VarArg;
+       GmfMshSct *msh;
+       KwdSct *kwd;
+
+       if( (MshIdx < 1) || (MshIdx > MaxMsh) )
+               return(0);
+
+       msh = GmfMshTab[ MshIdx ];
+       RecBlk(msh, msh->buf, 0);
+
+       if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) )
+               return(0);
+
+       kwd = &msh->KwdTab[ KwdCod ];
+
+       /* Read further arguments if this kw has a header */
+
+       if(strlen(GmfKwdFmt[ KwdCod ][2]))
+       {
+               va_start(VarArg, KwdCod);
+               NmbLin = va_arg(VarArg, int);
+
+               if(!strcmp(GmfKwdFmt[ KwdCod ][3], "sr"))
+               {
+                       kwd->NmbTyp = va_arg(VarArg, int);
+                       TypTab = va_arg(VarArg, int *);
+
+                       for(i=0;i<kwd->NmbTyp;i++)
+                               kwd->TypTab[i] = TypTab[i];
+               }
+
+               va_end(VarArg);
+       }
+
+       /* Setup the kwd info */
+
+       ExpFmt(msh, KwdCod);
+
+       if(!kwd->typ)
+               return(0);
+       else if(kwd->typ == InfKwd)
+               kwd->NmbLin = 1;
+       else
+               kwd->NmbLin = NmbLin;
+
+       /* Store the next kwd position in binary file */
+
+       if( (msh->typ & Bin) && msh->NexKwdPos )
+       {
+               CurPos = ftell(msh->hdl);
+               fseek(msh->hdl, msh->NexKwdPos, SEEK_SET);
+               SetPos(msh, CurPos);
+               fseek(msh->hdl, CurPos, SEEK_SET);
+       }
+
+       /* Write the header */
+
+       if(msh->typ & Asc)
+       {
+               fprintf(msh->hdl, "\n%s\n", GmfKwdFmt[ KwdCod ][0]);
+
+               if(kwd->typ != InfKwd)
+                       fprintf(msh->hdl, "%d\n", kwd->NmbLin);
+
+               /* In case of solution field, write the extended header */
+
+               if(kwd->typ == SolKwd)
+               {
+                       fprintf(msh->hdl, "%d ", kwd->NmbTyp);
+
+                       for(i=0;i<kwd->NmbTyp;i++)
+                               fprintf(msh->hdl, "%d ", kwd->TypTab[i]);
+
+                       fprintf(msh->hdl, "\n\n");
+               }
+       }
+       else
+       {
+               RecWrd(msh, (unsigned char *)&KwdCod);
+               msh->NexKwdPos = ftell(msh->hdl);
+               SetPos(msh, 0);
+
+               if(kwd->typ != InfKwd)
+                       RecWrd(msh, (unsigned char *)&kwd->NmbLin);
+
+               /* In case of solution field, write the extended header at once */
+
+               if(kwd->typ == SolKwd)
+               {
+                       RecWrd(msh, (unsigned char *)&kwd->NmbTyp);
+
+                       for(i=0;i<kwd->NmbTyp;i++)
+                               RecWrd(msh, (unsigned char *)&kwd->TypTab[i]);
+               }
+       }
+       msh->pos = 0;
+       return(kwd->NmbLin);
+}
+
+
+/*----------------------------------------------------------*/
+/* Read a full line from the current kwd                                       */
+/*----------------------------------------------------------*/
+
+void GmfGetLin(int MshIdx, int KwdCod, ...)
+{
+       int i, j;
+       float *FltSolTab;
+       double *DblSolTab;
+       va_list VarArg;
+       GmfMshSct *msh = GmfMshTab[ MshIdx ];
+       KwdSct *kwd = &msh->KwdTab[ KwdCod ];
+
+       /* Start decoding the arguments */
+
+       va_start(VarArg, KwdCod);
+
+       if(kwd->typ != SolKwd)
+       {
+               if(msh->ver == 1)
+               {
+                       if(msh->typ & Asc)
+                       {
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               fscanf(msh->hdl, "%f", va_arg(VarArg, float *));
+                                       else
+                                               fscanf(msh->hdl, "%d", va_arg(VarArg, int *));
+                       }
+                       else
+                       {
+                               ScaBlk(msh, msh->buf, kwd->SolSiz);
+
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               *(va_arg(VarArg, float *)) = msh->FltBuf[i];
+                                       else
+                                               *(va_arg(VarArg, int *)) = msh->IntBuf[i];
+                       }
+               }
+               else
+               {
+                       if(msh->typ & Asc)
+                       {
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               fscanf(msh->hdl, "%lf", va_arg(VarArg, double *));
+                                       else
+                                               fscanf(msh->hdl, "%d", va_arg(VarArg, int *));
+                       }
+                       else
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               ScaDblWrd(msh, (unsigned char *)va_arg(VarArg, double *));
+                                       else
+                                               ScaWrd(msh, (unsigned char *)va_arg(VarArg, int *));
+               }
+       }
+       else
+       {
+               if(msh->ver == 1)
+               {
+                       FltSolTab = va_arg(VarArg, float *);
+
+                       if(msh->typ & Asc)
+                               for(j=0;j<kwd->SolSiz;j++)
+                                       fscanf(msh->hdl, "%f", &FltSolTab[j]);
+                       else
+                               ScaBlk(msh, (unsigned char *)FltSolTab, kwd->NmbWrd);
+               }
+               else
+               {
+                       DblSolTab = va_arg(VarArg, double *);
+
+                       if(msh->typ & Asc)
+                               for(j=0;j<kwd->SolSiz;j++)
+                                       fscanf(msh->hdl, "%lf", &DblSolTab[j]);
+                       else
+                               for(j=0;j<kwd->SolSiz;j++)
+                                       ScaDblWrd(msh, (unsigned char *)&DblSolTab[j]);
+               }
+       }
+
+       va_end(VarArg);
+}
+
+
+/*----------------------------------------------------------*/
+/* Write a full line from the current kwd                                      */
+/*----------------------------------------------------------*/
+
+void GmfSetLin(int MshIdx, int KwdCod, ...)
+{
+       int i, j, pos, *IntBuf;
+       float *FltSolTab;
+       double *DblSolTab, *DblBuf;
+       va_list VarArg;
+       GmfMshSct *msh = GmfMshTab[ MshIdx ];
+       KwdSct *kwd = &msh->KwdTab[ KwdCod ];
+
+       /* Start decoding the arguments */
+
+       va_start(VarArg, KwdCod);
+
+       if(kwd->typ != SolKwd)
+       {
+               if(msh->ver == 1)
+               {
+                       if(msh->typ & Asc)
+                       {
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               fprintf(msh->hdl, "%g ", (float)va_arg(VarArg, double));
+                                       else
+                                               fprintf(msh->hdl, "%d ", va_arg(VarArg, int));
+                       }
+                       else
+                       {
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               msh->FltBuf[i] = va_arg(VarArg, double);
+                                       else
+                                               msh->IntBuf[i] = va_arg(VarArg, int);
+
+                               RecBlk(msh, msh->buf, kwd->SolSiz);
+                       }
+               }
+               else
+               {
+                       if(msh->typ & Asc)
+                       {
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                               fprintf(msh->hdl, "%.15lg ", va_arg(VarArg, double));
+                                       else
+                                               fprintf(msh->hdl, "%d ", va_arg(VarArg, int));
+                       }
+                       else
+                       {
+                               pos = 0;
+
+                               for(i=0;i<kwd->SolSiz;i++)
+                                       if(kwd->fmt[i] == 'r')
+                                       {
+                                               DblBuf = (void *)&msh->buf[ pos ];
+                                               *DblBuf = va_arg(VarArg, double);
+                                               pos += 8;
+                                       }
+                                       else
+                                       {
+                                               IntBuf = (void *)&msh->buf[ pos ];
+                                               *IntBuf = va_arg(VarArg, int);
+                                               pos += 4;
+                                       }
+                               RecBlk(msh, msh->buf, kwd->NmbWrd);
+                       }
+               }
+       }
+       else
+       {
+               if(msh->ver == 1)
+               {
+                       FltSolTab = va_arg(VarArg, float *);
+
+                       if(msh->typ & Asc)
+                               for(j=0;j<kwd->SolSiz;j++)
+                                       fprintf(msh->hdl, "%g ", FltSolTab[j]);
+                       else
+                               RecBlk(msh, (unsigned char *)FltSolTab, kwd->NmbWrd);
+               }
+               else
+               {
+                       DblSolTab = va_arg(VarArg, double *);
+
+                       if(msh->typ & Asc)
+                               for(j=0;j<kwd->SolSiz;j++)
+                                       fprintf(msh->hdl, "%.15lg ", DblSolTab[j]);
+                       else
+                               RecBlk(msh, (unsigned char *)DblSolTab, kwd->NmbWrd);
+               }
+       }
+
+       va_end(VarArg);
+
+       if(msh->typ & Asc)
+               fprintf(msh->hdl, "\n");
+}
+
+
+/*----------------------------------------------------------*/
+/* Private procedure for transmesh : copy a whole line         */
+/*----------------------------------------------------------*/
+
+void GmfCpyLin(int InpIdx, int OutIdx, int KwdCod)
+{
+       double d;
+       float f;
+       int i, a;
+       GmfMshSct *InpMsh = GmfMshTab[ InpIdx ], *OutMsh = GmfMshTab[ OutIdx ];
+       KwdSct *kwd = &InpMsh->KwdTab[ KwdCod ];
+
+       for(i=0;i<kwd->SolSiz;i++)
+       {
+               if(kwd->fmt[i] == 'r')
+               {
+                       if(InpMsh->ver == 1)
+                       {
+                               if(InpMsh->typ & Asc)
+                                       fscanf(InpMsh->hdl, "%f", &f);
+                               else
+                                       ScaWrd(InpMsh, (unsigned char *)&f);
+
+                               d = f;
+                       }
+                       else
+                       {
+                               if(InpMsh->typ & Asc)
+                                       fscanf(InpMsh->hdl, "%lf", &d);
+                               else
+                                       ScaDblWrd(InpMsh, (unsigned char *)&d);
+
+                               f = (float)d;
+                       }
+
+                       if(OutMsh->ver == 1)
+                               if(OutMsh->typ & Asc)
+                                       fprintf(OutMsh->hdl, "%g ", f);
+                               else
+                                       RecWrd(OutMsh, (unsigned char *)&f);
+                       else
+                               if(OutMsh->typ & Asc)
+                                       fprintf(OutMsh->hdl, "%.15g ", d);
+                               else
+                                       RecDblWrd(OutMsh, (unsigned char *)&d);
+               }
+               else
+               {
+                       if(InpMsh->typ & Asc)
+                               fscanf(InpMsh->hdl, "%d", &a);
+                       else
+                               ScaWrd(InpMsh, (unsigned char *)&a);
+
+                       if(OutMsh->typ & Asc)
+                               fprintf(OutMsh->hdl, "%d ", a);
+                       else
+                               RecWrd(OutMsh, (unsigned char *)&a);
+               }
+       }
+
+       if(OutMsh->typ & Asc)
+               fprintf(OutMsh->hdl, "\n");
+}
+
+
+/*----------------------------------------------------------*/
+/* Find every kw present in a meshfile                                         */
+/*----------------------------------------------------------*/
+
+static int ScaKwdTab(GmfMshSct *msh)
+{
+       int KwdCod;
+       long  NexPos, CurPos, EndPos;
+       char str[ GmfStrSiz ];
+
+       if(msh->typ & Asc)
+       {
+               /* Scan each string in the file until the end */
+
+               while(fscanf(msh->hdl, "%s", str) != EOF)
+               {
+                       /* Fast test in order to reject quickly the numeric values */
+
+                       if(isalpha(str[0]))
+                       {
+                               /* Search which kwd code this string is associated with, 
+                                       then get its header and save the curent position in file (just before the data) */
+
+                               for(KwdCod=1; KwdCod<= GmfMaxKwd; KwdCod++)
+                                       if(!strcmp(str, GmfKwdFmt[ KwdCod ][0]))
+                                       {
+                                               ScaKwdHdr(msh, KwdCod);
+                                               break;
+                                       }
+                       }
+                       else if(str[0] == '#')
+                               while(fgetc(msh->hdl) != '\n');
+               }
+       }
+       else
+       {
+               /* Get file size */
+
+               CurPos = ftell(msh->hdl);
+               fseek(msh->hdl, 0, SEEK_END);
+               EndPos = ftell(msh->hdl);
+               fseek(msh->hdl, CurPos, SEEK_SET);
+
+               /* Jump through kwd positions in the file */
+
+               do
+               {
+                       /* Get the kwd code and the next kwd position */
+
+                       ScaWrd(msh, (unsigned char *)&KwdCod);
+                       NexPos = GetPos(msh);
+
+                       if(NexPos > EndPos)
+                               return(0);
+
+                       /* Check if this kwd belongs to this mesh version */
+
+                       if( (KwdCod >= 1) && (KwdCod <= GmfMaxKwd) )
+                               ScaKwdHdr(msh, KwdCod);
+
+                       /* Go to the next kwd */
+
+                       if(NexPos)
+                               fseek(msh->hdl, NexPos, SEEK_SET);
+               }while(NexPos && (KwdCod != GmfEnd));
+       }
+
+       return(1);
+}
+
+
+/*----------------------------------------------------------*/
+/* Read and setup the keyword's header                                         */
+/*----------------------------------------------------------*/
+
+static void ScaKwdHdr(GmfMshSct *msh, int KwdCod)
+{
+       int i;
+       KwdSct *kwd = &msh->KwdTab[ KwdCod ];
+
+       if(!strcmp("i", GmfKwdFmt[ KwdCod ][2]))
+       {
+               if(msh->typ & Asc)
+                       fscanf(msh->hdl, "%d", &kwd->NmbLin);
+               else
+                       ScaWrd(msh, (unsigned char *)&kwd->NmbLin);
+       }
+       else
+               kwd->NmbLin = 1;
+
+       if(!strcmp("sr", GmfKwdFmt[ KwdCod ][3]))
+       {
+               if(msh->typ & Asc)
+               {
+                       fscanf(msh->hdl, "%d", &kwd->NmbTyp);
+
+                       for(i=0;i<kwd->NmbTyp;i++)
+                               fscanf(msh->hdl, "%d", &kwd->TypTab[i]);
+               }
+               else
+               {
+                       ScaWrd(msh, (unsigned char *)&kwd->NmbTyp);
+
+                       for(i=0;i<kwd->NmbTyp;i++)
+                               ScaWrd(msh, (unsigned char *)&kwd->TypTab[i]);
+               }
+       }
+
+       ExpFmt(msh, KwdCod);
+       kwd->pos = ftell(msh->hdl);
+}
+
+
+/*----------------------------------------------------------*/
+/* Expand the compacted format and compute the line size       */
+/*----------------------------------------------------------*/
+
+static void ExpFmt(GmfMshSct *msh, int KwdCod)
+{
+       int i, j, TmpSiz=0;
+       char chr, *InpFmt = GmfKwdFmt[ KwdCod ][3];
+       KwdSct *kwd = &msh->KwdTab[ KwdCod ];
+
+       /* Set the kwd's type */
+
+       if(!strlen(GmfKwdFmt[ KwdCod ][2]))
+               kwd->typ = InfKwd;
+       else if(!strcmp(InpFmt, "sr"))
+               kwd->typ = SolKwd;
+       else
+               kwd->typ = RegKwd;
+
+       /* Get the solution-field's size */
+
+       if(kwd->typ == SolKwd)
+               for(i=0;i<kwd->NmbTyp;i++)
+                       switch(kwd->TypTab[i])
+                       {
+                               case GmfSca    : TmpSiz += 1; break;
+                               case GmfVec    : TmpSiz += msh->dim; break;
+                               case GmfSymMat : TmpSiz += (msh->dim * (msh->dim+1)) / 2; break;
+                               case GmfMat    : TmpSiz += msh->dim * msh->dim; break;
+                       }
+
+       /* Scan each character from the format string */
+
+       i = kwd->SolSiz = kwd->NmbWrd = 0;
+
+       while(i < strlen(InpFmt))
+       {
+               chr = InpFmt[ i++ ];
+
+               if(chr == 'd')
+               {
+                       chr = InpFmt[i++];
+
+                       for(j=0;j<msh->dim;j++)
+                               kwd->fmt[ kwd->SolSiz++ ] = chr;
+               }
+               else if(chr == 's')
+               {
+                       chr = InpFmt[i++];
+
+                       for(j=0;j<TmpSiz;j++)
+                               kwd->fmt[ kwd->SolSiz++ ] = chr;
+               }
+               else
+                       kwd->fmt[ kwd->SolSiz++ ] = chr;
+       }
+
+       for(i=0;i<kwd->SolSiz;i++)
+               if(kwd->fmt[i] == 'i')
+                       kwd->NmbWrd++;
+               else if(msh->ver >= 2)
+                       kwd->NmbWrd += 2;
+               else
+                       kwd->NmbWrd++;
+}
+
+
+/*----------------------------------------------------------*/
+/* Read a four bytes word from a mesh file                                     */
+/*----------------------------------------------------------*/
+
+static void ScaWrd(GmfMshSct *msh, unsigned char *wrd)
+{
+       unsigned char swp;
+
+       fread(wrd, WrdSiz, 1, msh->hdl);
+
+       if(msh->cod == 1)
+               return;
+
+       swp = wrd[3];
+       wrd[3] = wrd[0];
+       wrd[0] = swp;
+
+       swp = wrd[2];
+       wrd[2] = wrd[1];
+       wrd[1] = swp;
+}
+
+
+/*----------------------------------------------------------*/
+/* Read an eight bytes word from a mesh file                           */
+/*----------------------------------------------------------*/
+
+static void ScaDblWrd(GmfMshSct *msh, unsigned char *wrd)
+{
+       int i;
+       unsigned char swp;
+
+       fread(wrd, WrdSiz, 2, msh->hdl);
+
+       if(msh->cod == 1)
+               return;
+
+       for(i=0;i<4;i++)
+       {
+               swp = wrd[7-i];
+               wrd[7-i] = wrd[i];
+               wrd[i] = swp;
+       }
+}
+
+
+/*----------------------------------------------------------*/
+/* Read ablock of four bytes word from a mesh file                     */
+/*----------------------------------------------------------*/
+
+static void ScaBlk(GmfMshSct *msh, unsigned char *blk, int siz)
+{
+       int i, j;
+       unsigned char swp, *wrd;
+
+       fread(blk, WrdSiz, siz, msh->hdl);
+
+       if(msh->cod == 1)
+               return;
+
+       for(i=0;i<siz;i++)
+       {
+               wrd = &blk[ i * 4 ];
+
+               for(j=0;j<2;j++)
+               {
+                       swp = wrd[ 3-j ];
+                       wrd[ 3-j ] = wrd[j];
+                       wrd[j] = swp;
+               }
+       }
+}
+
+
+/*----------------------------------------------------------*/
+/* Read a 4 or 8 bytes position in mesh file                           */
+/*----------------------------------------------------------*/
+
+static long GetPos(GmfMshSct *msh)
+{
+       int IntVal;
+       long pos;
+
+       if(msh->ver >= 3)
+               ScaDblWrd(msh, (unsigned char*)&pos);
+       else
+       {
+               ScaWrd(msh, (unsigned char*)&IntVal);
+               pos = IntVal;
+       }
+
+       return(pos);
+}
+
+
+/*----------------------------------------------------------*/
+/* Write a four bytes word to a mesh file                                      */
+/*----------------------------------------------------------*/
+
+static void RecWrd(GmfMshSct *msh, unsigned char *wrd)
+{
+       fwrite(wrd, WrdSiz, 1, msh->hdl);
+}
+
+
+/*----------------------------------------------------------*/
+/* Write an eight bytes word to a mesh file                                    */
+/*----------------------------------------------------------*/
+
+static void RecDblWrd(GmfMshSct *msh, unsigned char *wrd)
+{
+       fwrite(wrd, WrdSiz, 2, msh->hdl);
+}
+
+
+/*----------------------------------------------------------*/
+/* Write a block of four bytes word to a mesh file                     */
+/*----------------------------------------------------------*/
+
+static void RecBlk(GmfMshSct *msh, unsigned char *blk, int siz)
+{
+       /* Copy this line-block into the main mesh buffer */
+
+       if(siz)
+       {
+               memcpy(&msh->blk[ msh->pos ], blk, siz * WrdSiz);
+               msh->pos += siz * WrdSiz;
+       }
+
+       /* When the buffer is full or this procedure is called with a 0 size, flush the cache on disk */
+
+       if( (msh->pos > BufSiz) || (!siz && msh->pos) )
+       {
+               fwrite(msh->blk, 1, msh->pos, msh->hdl);
+               msh->pos = 0;
+       }
+}
+
+
+/*----------------------------------------------------------*/
+/* Read a 4 or 8 bytes position in mesh file                           */
+/*----------------------------------------------------------*/
+
+static void SetPos(GmfMshSct *msh, long pos)
+{
+       int IntVal;
+
+       if(msh->ver >= 3)
+               RecDblWrd(msh, (unsigned char*)&pos);
+       else
+       {
+               IntVal = pos;
+               RecWrd(msh, (unsigned char*)&IntVal);
+       }
+}
diff --git a/src/GHS3DPlugin/libmesh5.h b/src/GHS3DPlugin/libmesh5.h
new file mode 100755 (executable)
index 0000000..cfaf9fb
--- /dev/null
@@ -0,0 +1,152 @@
+
+
+/*----------------------------------------------------------*/
+/*                                                                                                                     */
+/*                                             LIBMESH V 5.45                                          */
+/*                                                                                                                     */
+/*----------------------------------------------------------*/
+/*                                                                                                                     */
+/*     Description:            handle .meshb file format I/O           */
+/*     Author:                         Loic MARECHAL                                           */
+/*     Creation date:          feb 16 2007                                                     */
+/*     Last modification:      sep 27 2010                                                     */
+/*                                                                                                                     */
+/*----------------------------------------------------------*/
+
+
+/*----------------------------------------------------------*/
+/* Defines                                                                                                     */
+/*----------------------------------------------------------*/
+
+#define GmfStrSiz 1024
+#define GmfMaxTyp 1000
+#define GmfMaxKwd 79
+#define GmfMshVer 1
+#define GmfRead 1
+#define GmfWrite 2
+#define GmfSca 1
+#define GmfVec 2
+#define GmfSymMat 3
+#define GmfMat 4
+#define GmfFloat 1
+#define GmfDouble 2
+
+enum GmfKwdCod
+{
+       GmfReserved1, \
+       GmfVersionFormatted, \
+       GmfReserved2, \
+       GmfDimension, \
+       GmfVertices, \
+       GmfEdges, \
+       GmfTriangles, \
+       GmfQuadrilaterals, \
+       GmfTetrahedra, \
+       GmfPrisms, \
+       GmfHexahedra, \
+       GmfIterationsAll, \
+       GmfTimesAll, \
+       GmfCorners, \
+       GmfRidges, \
+       GmfRequiredVertices, \
+       GmfRequiredEdges, \
+       GmfRequiredTriangles, \
+       GmfRequiredQuadrilaterals, \
+       GmfTangentAtEdgeVertices, \
+       GmfNormalAtVertices, \
+       GmfNormalAtTriangleVertices, \
+       GmfNormalAtQuadrilateralVertices, \
+       GmfAngleOfCornerBound, \
+       GmfTrianglesP2, \
+       GmfEdgesP2, \
+       GmfSolAtPyramids, \
+       GmfQuadrilateralsQ2, \
+       GmfISolAtPyramids, \
+       GmfReserved6, \
+       GmfTetrahedraP2, \
+       GmfReserved7, \
+       GmfReserved8, \
+       GmfHexahedraQ2, \
+       GmfReserved9, \
+       GmfReserved10, \
+       GmfReserved17, \
+       GmfReserved18, \
+       GmfReserved19, \
+       GmfReserved20, \
+       GmfReserved21, \
+       GmfReserved22, \
+       GmfReserved23, \
+       GmfReserved24, \
+       GmfReserved25, \
+       GmfReserved26, \
+       GmfPolyhedra, \
+       GmfPolygons, \
+       GmfReserved29, \
+       GmfPyramids, \
+       GmfBoundingBox, \
+       GmfBody, \
+       GmfPrivateTable, \
+       GmfReserved33, \
+       GmfEnd, \
+       GmfReserved34, \
+       GmfReserved35, \
+       GmfReserved36, \
+       GmfReserved37, \
+       GmfTangents, \
+       GmfNormals, \
+       GmfTangentAtVertices, \
+       GmfSolAtVertices, \
+       GmfSolAtEdges, \
+       GmfSolAtTriangles, \
+       GmfSolAtQuadrilaterals, \
+       GmfSolAtTetrahedra, \
+       GmfSolAtPrisms, \
+       GmfSolAtHexahedra, \
+       GmfDSolAtVertices, \
+       GmfISolAtVertices, \
+       GmfISolAtEdges, \
+       GmfISolAtTriangles, \
+       GmfISolAtQuadrilaterals, \
+       GmfISolAtTetrahedra, \
+       GmfISolAtPrisms, \
+       GmfISolAtHexahedra, \
+       GmfIterations, \
+       GmfTime, \
+       GmfReserved38
+};
+
+
+/*----------------------------------------------------------*/
+/* External procedures                                                                         */
+/*----------------------------------------------------------*/
+
+extern int GmfOpenMesh(const char *, int, ...);
+extern int GmfCloseMesh(int);
+extern int GmfStatKwd(int, int, ...);
+extern int GmfGotoKwd(int, int);
+extern int GmfSetKwd(int, int, ...);
+extern void GmfGetLin(int, int, ...);
+extern void GmfSetLin(int, int, ...);
+
+
+/*----------------------------------------------------------*/
+/* Fortran 77 API                                                                                      */
+/*----------------------------------------------------------*/
+
+#if defined(F77_NO_UNDER_SCORE)
+#define call(x) x
+#else
+#define call(x) x ## _
+#endif
+
+
+/*----------------------------------------------------------*/
+/* Transmesh private API                                                                       */
+/*----------------------------------------------------------*/
+
+#ifdef TRANSMESH
+
+extern char *GmfKwdFmt[ GmfMaxKwd + 1 ][4];
+extern int GmfCpyLin(int, int, int);
+
+#endif