]> SALOME platform Git repositories - plugins/blsurfplugin.git/commitdiff
Salome HOME
Merge from V6_6_BR (V6_6_0rc2) 11/12/2012
authorvsr <vsr@opencascade.com>
Tue, 11 Dec 2012 12:03:18 +0000 (12:03 +0000)
committervsr <vsr@opencascade.com>
Tue, 11 Dec 2012 12:03:18 +0000 (12:03 +0000)
adm_local/cmake_files/FindCADSURF.cmake
configure.ac
doc/salome/gui/BLSURFPLUGIN/input/blsurf_hypo.doc
resources/BLSURFPlugin.xml
src/BLSURFPlugin/BLSURFPlugin_BLSURF.cxx
src/BLSURFPlugin/BLSURFPlugin_BLSURF.hxx
src/GUI/BLSURFPluginGUI_Dlg.h
src/GUI/Makefile.am

index 1954cb732e610e22ae4544199068bd6c0b3b05fd..6822a40296b9b7b2b6a2d65439f6226e1a077959 100644 (file)
@@ -19,8 +19,8 @@
 
 SET(MESHGEMSHOME $ENV{MESHGEMSHOME})
 FIND_PATH(CADSURF_INCLUDES_DIR meshgems/cadsurf.h ${MESHGEMSHOME}/include)
-SET(CADSURF_INCLUDES)
-SET(CADSURF_INCLUDES ${CADSURF_INCLUDES} -I${CADSURF_INCLUDES_DIR})
+SET(MESHGEMS_CADSURF_INCLUDES)
+SET(MESHGEMS_CADSURF_INCLUDES ${MESHGEMS_CADSURF_INCLUDES} -I${CADSURF_INCLUDES_DIR})
 
 SET(CADSURF_LIBS_PATHS)
 SET(CADSURF_LIBS_PATHS ${CADSURF_LIBS_PATHS} ${MESHGEMSHOME}/lib)
@@ -38,9 +38,9 @@ FIND_LIBRARY(mg-cadsurf mg-cadsurf PATHS ${CADSURF_LIBS_PATHS})
 FIND_LIBRARY(mg-precad mg-precad PATHS ${CADSURF_LIBS_PATHS})
 FIND_LIBRARY(meshgems meshgems PATHS ${CADSURF_LIBS_PATHS})
 
-SET(CADSURF_LIBS)
-SET(CADSURF_LIBS ${CADSURF_LIBS} ${mg-cadsurf})
-SET(CADSURF_LIBS ${CADSURF_LIBS} ${mg-precad})
+SET(MESHGEMS_CADSURF_LIBS)
+SET(MESHGEMS_CADSURF_LIBS ${MESHGEMS_CADSURF_LIBS} ${mg-cadsurf})
+SET(MESHGEMS_CADSURF_LIBS ${MESHGEMS_CADSURF_LIBS} ${mg-precad})
 IF(meshgems)
-  SET(CADSURF_LIBS ${CADSURF_LIBS} ${meshgems})
+  SET(MESHGEMS_CADSURF_LIBS ${MESHGEMS_CADSURF_LIBS} ${meshgems})
 ENDIF(meshgems)
index d68d4fd523cc0dda7ccebd2f1a0d8f031799ea54..3580666470b34875f0a4a4365ec7e4de63fefb33 100644 (file)
@@ -23,7 +23,7 @@
 # Author : Vadim SANDLER, Open CASCADE S.A.S (vadim.sandler@opencascade.com)
 # ---
 #
-AC_INIT([Salome2 Project BLSURFPLUGIN module], [6.5.0], [webmaster.salome@opencascade.com], [SalomeBLSURFPLUGIN])
+AC_INIT([Salome2 Project BLSURFPLUGIN module], [6.6.0], [webmaster.salome@opencascade.com], [SalomeBLSURFPLUGIN])
 AC_CONFIG_AUX_DIR(adm_local/unix/config_files)
 AC_CANONICAL_HOST
 AC_CANONICAL_TARGET
index 951745a691769425074bbc8024e26eb2a8635fda..481510c14696430db33aa7bff661aca9768a0f88 100644 (file)
@@ -54,12 +54,12 @@ The smaller this distance is, the closer the mesh is to the exact surface (only
 - <b>Anisotropic</b> - if checked, this parameter defines the maximum anisotropic ratio of the metric governing the anisotropic meshing process.
 The default value (0) means that the metric (and thus the generated elements) can be arbitrarily stretched.
 
-- <b>Remove tiny edges</b> - if checked, the bad elements (slivers) are removed from the generated mesh.
-The bad element value defines the aspect ratio triggering the "bad element” classification.
-
-- <b>Remove bad elements</b> - if checked, the tiny (nano) edges are removed from the generated mesh.
+- <b>Remove tiny edges</b> - if checked, the tiny (nano) edges are removed from the generated mesh.
 The tiny edge value defines the minimal length under which an edge is considered to be a tiny one.
 
+- <b>Remove bad elements</b> - if checked, the bad elements (slivers) are removed from the generated mesh.
+The bad element value defines the aspect ratio triggering the "bad element” classification.
+
 - <b>Mesh optimisation</b> - if checked, the mesh will be optimized in order to get better shaped elements.
 
 - <b>Allow Quadrangles</b> - if checked, allows the creation of quadrilateral elements.
index 857d339f5fec98382515613372062404a2394c73..93cbaf79fc96a1f69ef059b125eef734b84c4a7f 100644 (file)
@@ -73,6 +73,7 @@
       <python-wrap>
         <algo>BLSURF=Triangle(algo=smesh.BLSURF)</algo>
         <hypo>BLSURF_Parameters=Parameters()</hypo>
+        <hypo>ViscousLayers2D=ViscousLayers2D(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetIgnoreEdges())</hypo>
       </python-wrap>
     </algorithm>
 
index 173a7e1df3a6d42abecb883ba3696b0de3f72f03..01cc508e2323f2e8215d88c1dd5fcce80e9d0a5a 100644 (file)
@@ -58,47 +58,42 @@ extern "C"{
 #include <cstdlib>
 
 // OPENCASCADE includes
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakePolygon.hxx>
+//#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <BRepGProp.hxx>
+#include <BRepTools.hxx>
+#include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
+#include <GProp_GProps.hxx>
+#include <Geom2d_Curve.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Surface.hxx>
+#include <NCollection_Map.hxx>
+#include <Standard_ErrorHandler.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
+#include <TopTools_DataMapOfShapeInteger.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
-#include <NCollection_Map.hxx>
-
-#include <Geom_Surface.hxx>
-#include <Handle_Geom_Surface.hxx>
-#include <Geom2d_Curve.hxx>
-#include <Handle_Geom2d_Curve.hxx>
-#include <Geom_Curve.hxx>
-#include <Handle_Geom_Curve.hxx>
-#include <Handle_AIS_InteractiveObject.hxx>
-#include <TopoDS_Vertex.hxx>
+#include <TopoDS_Compound.hxx>
 #include <TopoDS_Edge.hxx>
-#include <TopoDS_Wire.hxx>
 #include <TopoDS_Face.hxx>
-
-#include <gp_Pnt2d.hxx>
-#include <gp_Pnt.hxx>
-#include <TopTools_IndexedMapOfShape.hxx>
 #include <TopoDS_Shape.hxx>
-#include <BRep_Builder.hxx>
-#include <BRepBuilderAPI_MakeVertex.hxx>
-#include <BRepTools.hxx>
-
-#include <TopTools_DataMapOfShapeInteger.hxx>
-#include <GProp_GProps.hxx>
-#include <BRepGProp.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <TopoDS_Wire.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Pnt2d.hxx>
+#include <gp_XY.hxx>
+#include <gp_XYZ.hxx>
 
 #ifndef WNT
 #include <fenv.h>
 #endif
 
-#include <Standard_ErrorHandler.hxx>
-#include <GeomAPI_ProjectPointOnCurve.hxx>
-#include <GeomAPI_ProjectPointOnSurf.hxx>
-#include <gp_XY.hxx>
-#include <gp_XYZ.hxx>
-#include <TopTools_MapOfShape.hxx>
-
 /* ==================================
  * ===========  PYTHON ==============
  * ==================================*/
@@ -1246,83 +1241,12 @@ namespace
   };
 
   // --------------------------------------------------------------------------
-  /*!
-   * \brief A curve of internal boundary of viscous layer
-   */
-  class TProxyPCurve : public Geom2d_Curve
+  // comparator to sort nodes and sub-meshes
+  struct ShapeTypeCompare
   {
-    const UVPtStructVec& _noDataVec;
-    Handle(Geom2d_Curve) _pcurve;
-    double               _sign;
-
-    int locateU( const double U, double& r ) const
-    {
-      const double eps = 1e-10;
-      double normU =
-        ( U - _noDataVec[0].param ) / ( _noDataVec.back().param - _noDataVec[0].param );
-      int i = int( normU * ( _noDataVec.size() - 1 ));
-      while ( _noDataVec[ i ].normParam > normU + eps && i != 0 )
-        --i;
-      while ( i+1 < _noDataVec.size() && _noDataVec[ i+1 ].normParam < normU - eps )
-        ++i;
-      //cout << "U " << U << " normU " << normU << " i " << i << " _noDataVec.size() " << _noDataVec.size() << endl;
-      r = 0;
-      if ( normU > eps && normU < 1.- eps && i+1 < _noDataVec.size() )
-        r = ( normU - _noDataVec[i].normParam ) /
-          ( _noDataVec[i+1].normParam - _noDataVec[i].normParam );
-      return i;
-    }
-  public:
-    TProxyPCurve( const SMESH_ProxyMesh::SubMesh* edgeSM,
-                  Handle(Geom2d_Curve)            pcurve)
-      : _noDataVec( edgeSM->GetUVPtStructVec() ), _pcurve( pcurve )
-    {
-      _sign = ( _noDataVec.back().param > _noDataVec[0].param ? +1 : -1 );
-    }
-    void D0(const Standard_Real U, gp_Pnt2d& P) const
-    {
-      double r;
-      int i = locateU( U, r );
-      P.SetCoord( _noDataVec[i].u, _noDataVec[i].v );
-      if ( r > 0 )
-        P = ( 1-r ) * P.XY() + r * gp_XY( _noDataVec[i+1].u, _noDataVec[i+1].v );
-      //cout << "U " << U << "i,r = ( "<< i << ", "<< r << " )\t P ( " << P.X() << ", " << P.Y() <<" )" << endl;
-    }
-    gp_Vec2d DN(const Standard_Real U,const Standard_Integer N) const
-    {
-      double r;
-      int i    = locateU( U, r );
-      double u = _noDataVec[i].param;
-      if ( r > 0 )
-        u = ( 1-r ) * u + r * _noDataVec[i+1].param;
-      gp_Vec2d dn = _pcurve->DN( u, N );
-      //cout << "U " << U << " D" << N << " " << dn.X() << ", " << dn.Y() << end;
-      return dn;
-    }
-    // virtual methods not used within BLSURF plug-in
-    void Reverse()                                         {}
-    Standard_Real ReversedParameter(const double U ) const { return 1 - U; }
-    Standard_Real FirstParameter() const                   { return 0;}
-    Standard_Real LastParameter() const                    { return 1;}
-    Standard_Boolean IsClosed() const                      { return false; }
-    Standard_Boolean IsPeriodic() const                    { return false; }
-    GeomAbs_Shape Continuity() const                       { return GeomAbs_G2; }
-    Standard_Boolean IsCN(const int N) const               { return false; }
-    void Transform(const gp_Trsf2d&)                       {}
-    Handle_Geom2d_Geometry Copy() const                    {return Handle_Geom2d_Geometry();}
-    void D1(const Standard_Real U,gp_Pnt2d& P,gp_Vec2d& V1) const {}
-    void D2(const Standard_Real U,gp_Pnt2d& P,gp_Vec2d& V1,gp_Vec2d& V2) const {}
-    void D3(const Standard_Real U,gp_Pnt2d& P,gp_Vec2d& V1,gp_Vec2d& V2,gp_Vec2d& V3) const {}
-
-    //DEFINE_STANDARD_RTTI(TProxyPCurve);
-  };
-  DEFINE_STANDARD_HANDLE(TProxyPCurve,Geom2d_Curve);
-  // --------------------------------------------------------------------------
-  // comparator to sort nodes by position in the following order:
-  // SMDS_TOP_FACE=2, SMDS_TOP_EDGE=1, SMDS_TOP_VERTEX=0, SMDS_TOP_3DSPACE=3
-  struct NodePosCompare
-  {
-    int operator()( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 )
+    // sort nodes by position in the following order:
+    // SMDS_TOP_FACE=2, SMDS_TOP_EDGE=1, SMDS_TOP_VERTEX=0, SMDS_TOP_3DSPACE=3
+    int operator()( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 ) const
     {
       SMDS_TypeOfPosition pos1 = n1->GetPosition()->GetTypeOfPosition();
       SMDS_TypeOfPosition pos2 = n2->GetPosition()->GetTypeOfPosition();
@@ -1330,6 +1254,15 @@ namespace
       if ( pos1 < pos2 || pos1 == SMDS_TOP_3DSPACE ) return 1;
       return -1;
     }
+    // sort sub-meshes in order: EDGE, VERTEX
+    bool operator()( const SMESHDS_SubMesh* s1, const SMESHDS_SubMesh* s2 ) const
+    {
+      int isVertex1 = ( s1 && s1->NbElements() == 0 );
+      int isVertex2 = ( s2 && s2->NbElements() == 0 );
+      if ( isVertex1 == isVertex2 )
+        return s1 < s2;
+      return isVertex1 < isVertex2;
+    }
   };
 
   //================================================================================
@@ -1382,7 +1315,7 @@ namespace
           }
           // make nodes created on the boundary of viscous layer replace nodes created
           // by BLSURF as their SMDS_Position is more correct
-          nodes.sort( NodePosCompare() );
+          nodes.sort( ShapeTypeCompare() );
           nodeGroupsToMerge.push_back( nodes );
         }
       }
@@ -1403,6 +1336,162 @@ namespace
     // cout << endl;
   }
 
+  //================================================================================
+  /*!
+   * \brief A temporary mesh used to compute mesh on a proxy FACE
+   */
+  //================================================================================
+
+  struct TmpMesh: public SMESH_Mesh
+  {
+    typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*, TIDCompare > TN2NMap;
+    TN2NMap     _tmp2origNN;
+    TopoDS_Face _proxyFace;
+
+    TmpMesh()
+    {
+      _myMeshDS = new SMESHDS_Mesh( _id, true );
+    }
+    //--------------------------------------------------------------------------------
+    /*!
+     * \brief Creates a FACE bound by viscous layers and mesh each its EDGE with 1 segment
+     */
+    //--------------------------------------------------------------------------------
+
+    const TopoDS_Face& makeProxyFace( SMESH_ProxyMesh::Ptr& viscousMesh,
+                                      const TopoDS_Face&    origFace)
+    {
+      // get data of nodes on inner boundary of viscous layers
+      SMESH_Mesh* origMesh = viscousMesh->GetMesh();
+      TError err;
+      TSideVector wireVec = StdMeshers_FaceSide::GetFaceWires(origFace, *origMesh,
+                                                              /*skipMediumNodes = */true,
+                                                              err, viscousMesh );
+      if ( err && err->IsKO() )
+        throw *err.get(); // it should be caught at SMESH_subMesh
+
+      // proxy nodes and corresponding tmp VERTEXes
+      std::vector<const SMDS_MeshNode*> origNodes;
+      std::vector<TopoDS_Vertex>        tmpVertex;
+
+      // create a proxy FACE
+      TopoDS_Shape origFaceCopy = origFace.EmptyCopied();
+      BRepBuilderAPI_MakeFace newFace( TopoDS::Face( origFaceCopy ));
+      for ( size_t iW = 0; iW != wireVec.size(); ++iW )
+      {
+        StdMeshers_FaceSidePtr& wireData = wireVec[iW];
+        const UVPtStructVec& wirePoints = wireData->GetUVPtStruct();
+        if ( wirePoints.size() < 3 )
+          continue;
+
+        BRepBuilderAPI_MakePolygon wire;
+        for ( size_t iN = 1; iN < wirePoints.size(); ++iN )
+        {
+          wire.Add( SMESH_TNodeXYZ( wirePoints[ iN ].node ));
+          origNodes.push_back( wirePoints[ iN ].node );
+          tmpVertex.push_back( wire.LastVertex() );
+        }
+        tmpVertex[0] = wire.FirstVertex();
+        wire.Close();
+        if ( !wire.IsDone() )
+          throw SALOME_Exception("BLSURFPlugin_BLSURF: BRepBuilderAPI_MakePolygon failed");
+        newFace.Add( wire );
+      }
+      _proxyFace = newFace;
+
+      // set a new shape to mesh
+      TopoDS_Compound auxCompoundToMesh;
+      BRep_Builder shapeBuilder;
+      shapeBuilder.MakeCompound( auxCompoundToMesh );
+      shapeBuilder.Add( auxCompoundToMesh, _proxyFace );
+      shapeBuilder.Add( auxCompoundToMesh, origMesh->GetShapeToMesh() );
+
+      ShapeToMesh( auxCompoundToMesh );
+
+      //TopExp_Explorer fExp( auxCompoundToMesh, TopAbs_FACE );
+      //_proxyFace = TopoDS::Face( fExp.Current() );
+
+
+      // Make input mesh for BLSURF: segments on EDGE's of newFace
+
+      // make nodes and fill in _tmp2origNN
+      //
+      SMESHDS_Mesh* tmpMeshDS = GetMeshDS();
+      for ( size_t i = 0; i < origNodes.size(); ++i )
+      {
+        GetSubMesh( tmpVertex[i] )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
+        if ( const SMDS_MeshNode* tmpN = SMESH_Algo::VertexNode( tmpVertex[i], tmpMeshDS ))
+          _tmp2origNN.insert( _tmp2origNN.end(), make_pair( tmpN, origNodes[i] ));
+        else
+          throw SALOME_Exception("BLSURFPlugin_BLSURF: a proxy vertex not meshed");
+      }
+
+      // make segments
+      TopoDS_Vertex v1, v2;
+      for ( TopExp_Explorer edge( _proxyFace, TopAbs_EDGE ); edge.More(); edge.Next() )
+      {
+        const TopoDS_Edge& E = TopoDS::Edge( edge.Current() );
+        TopExp::Vertices( E, v1, v2 );
+        const SMDS_MeshNode* n1 = SMESH_Algo::VertexNode( v1, tmpMeshDS );
+        const SMDS_MeshNode* n2 = SMESH_Algo::VertexNode( v2, tmpMeshDS );
+
+        if ( SMDS_MeshElement* seg = tmpMeshDS->AddEdge( n1, n2 ))
+          tmpMeshDS->SetMeshElementOnShape( seg, E );
+      }
+
+      return _proxyFace;
+    }
+
+    //--------------------------------------------------------------------------------
+    /*!
+     * \brief Fill in the origMesh with faces computed by BLSURF in this tmp mesh
+     */
+    //--------------------------------------------------------------------------------
+
+    void FillInOrigMesh( SMESH_Mesh&        origMesh,
+                         const TopoDS_Face& origFace )
+    {
+      SMESH_MesherHelper helper( origMesh );
+      helper.SetSubShape( origFace );
+      helper.SetElementsOnShape( true );
+
+      // iterate over tmp faces and copy them in origMesh
+      const SMDS_MeshNode* nodes[27];
+      const SMDS_MeshNode* nullNode = 0;
+      double xyz[3];
+      SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
+      while ( fIt->more() )
+      {
+        const SMDS_MeshElement* f = fIt->next();
+        SMDS_ElemIteratorPtr nIt = f->nodesIterator();
+        int nbN = 0;
+        for ( ; nIt->more(); ++nbN )
+        {
+          const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
+          TN2NMap::iterator n2nIt = 
+            _tmp2origNN.insert( _tmp2origNN.end(), make_pair( n, nullNode ));
+          if ( !n2nIt->second ) {
+            n->GetXYZ( xyz );
+            gp_XY uv = helper.GetNodeUV( _proxyFace, n );
+            n2nIt->second = helper.AddNode( xyz[0], xyz[1], xyz[2], uv.X(), uv.Y() );
+          }
+          nodes[ nbN ] = n2nIt->second;
+        }
+        switch( nbN ) {
+        case 3: helper.AddFace( nodes[0], nodes[1], nodes[2] ); break;
+        // case 6: helper.AddFace( nodes[0], nodes[1], nodes[2],
+        //                         nodes[3], nodes[4], nodes[5]); break;
+        case 4: helper.AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
+        // case 9: helper.AddFace( nodes[0], nodes[1], nodes[2], nodes[3],
+        //                         nodes[4], nodes[5], nodes[6], nodes[7], nodes[8]); break;
+        // case 8: helper.AddFace( nodes[0], nodes[1], nodes[2], nodes[3],
+        //                         nodes[4], nodes[5], nodes[6], nodes[7]); break;
+        }
+      }
+    }
+  };
+
+
 } // namespace
 
 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
@@ -1424,6 +1513,59 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
   // Fix problem with locales
   Kernel_Utils::Localizer aLocalizer;
 
+  if ( !compute( aMesh, aShape ))
+    return false;
+
+  if ( _haveViscousLayers )
+  {
+    // Compute viscous layers
+
+    TopTools_MapOfShape map;
+    for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
+    {
+      const TopoDS_Face& F = TopoDS::Face(face_iter.Current());
+      if ( !map.Add( F )) continue;
+      SMESH_ProxyMesh::Ptr viscousMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+      if ( !viscousMesh )
+        return false; // error in StdMeshers_ViscousLayers2D::Compute()
+
+      // Compute BLSURF mesh on viscous layers
+
+      if ( viscousMesh->NbProxySubMeshes() > 0 )
+      {
+        TmpMesh tmpMesh;
+        const TopoDS_Face& proxyFace = tmpMesh.makeProxyFace( viscousMesh, F );
+        if ( !compute( tmpMesh, proxyFace ))
+          return false;
+        tmpMesh.FillInOrigMesh( aMesh, F );
+      }
+    }
+
+    // Re-compute BLSURF mesh on the rest faces if the mesh was cleared
+
+    for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
+    {
+      const TopoDS_Face& F = TopoDS::Face(face_iter.Current());
+      SMESH_subMesh* fSM = aMesh.GetSubMesh( F );
+      if ( fSM->IsMeshComputed() ) continue;
+
+      if ( !compute( aMesh, aShape ))
+        return false;
+      break;
+    }
+  }
+  return true;
+}
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
+                                  const TopoDS_Shape& aShape)
+{
   /* create a distene context (generic object) */
   status_t status = STATUS_ERROR;
 
@@ -1434,10 +1576,9 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
   bool haveQuadraticSubMesh = SMESH_MesherHelper( aMesh ).IsQuadraticSubMesh( aShape );
   bool quadraticSubMeshAndViscousLayer = false;
   bool needMerge = false;
-  set< SMESHDS_SubMesh* > edgeSubmeshes, proxySubmeshes;
-  set< SMESHDS_SubMesh* >& mergeSubmeshes = edgeSubmeshes;
-
-  vector< SMESH_ProxyMesh::Ptr > viscousMeshPerFace;
+  typedef set< SMESHDS_SubMesh*, ShapeTypeCompare > TSubMeshSet;
+  TSubMeshSet edgeSubmeshes;
+  TSubMeshSet& mergeSubmeshes = edgeSubmeshes;
 
   TopTools_IndexedMapOfShape fmap;
   TopTools_IndexedMapOfShape emap;
@@ -1449,252 +1590,163 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
 #endif
 
-  for ( int iLoop = 0; iLoop < 2; ++iLoop ) /* mesh of the 1st loop will be
-                                               discarded if _haveViscousLayers */
-  {
-    context_t *ctx =  context_new();
+  context_t *ctx =  context_new();
 
-    /* Set the message callback in the working context */
-    context_set_message_callback(ctx, message_cb, &_comment);
+  /* Set the message callback in the working context */
+  context_set_message_callback(ctx, message_cb, &_comment);
 #ifdef WITH_SMESH_CANCEL_COMPUTE
-    _compute_canceled = false;
-    context_set_interrupt_callback(ctx, interrupt_cb, this);
+  _compute_canceled = false;
+  context_set_interrupt_callback(ctx, interrupt_cb, this);
 #else
-    context_set_interrupt_callback(ctx, interrupt_cb, NULL);
+  context_set_interrupt_callback(ctx, interrupt_cb, NULL);
 #endif
 
-    /* create the CAD object we will work on. It is associated to the context ctx. */
-    cad_t *c     = cad_new(ctx);
-    dcad_t *dcad = dcad_new(c);
-
-    FacesWithSizeMap.Clear();
-    FaceId2SizeMap.clear();
-    FaceId2ClassAttractor.clear();
-    FaceIndex2ClassAttractor.clear();
-    EdgesWithSizeMap.Clear();
-    EdgeId2SizeMap.clear();
-    VerticesWithSizeMap.Clear();
-    VertexId2SizeMap.clear();
-
-    /* Now fill the CAD object with data from your CAD
-     * environement. This is the most complex part of a successfull
-     * integration.
-     */
+  /* create the CAD object we will work on. It is associated to the context ctx. */
+  cad_t *c     = cad_new(ctx);
+  dcad_t *dcad = dcad_new(c);
 
-    // PreCAD
-    // If user requests it, send the CAD through Distene preprocessor : PreCAD
-    cad_t *cleanc = NULL; // preprocessed cad
-    precad_session_t *pcs = precad_session_new(ctx);
-    precad_data_set_cad(pcs, c);
-
-    cadsurf_session_t *css = cadsurf_session_new(ctx);
-
-    // an object that correctly deletes all cadsurf objects at destruction
-    BLSURF_Cleaner cleaner( ctx,css,c,dcad );
-
-    MESSAGE("BEGIN SetParameters");
-    bool use_precad = false;
-    SetParameters(
-                  // #if BLSURF_VERSION_LONG >= "3.1.1"
-                  //     c,
-                  // #endif
-                  _hypothesis, css, pcs, aMesh, &use_precad);
-    MESSAGE("END SetParameters");
-
-    haveQuadraticSubMesh = haveQuadraticSubMesh || (_hypothesis != NULL && _hypothesis->GetQuadraticMesh());
-    helper.SetIsQuadratic( haveQuadraticSubMesh );
-
-    // To remove as soon as quadratic mesh is allowed - BEGIN
-    // GDD: Viscous layer is not allowed with quadratic mesh
-    if (_haveViscousLayers && haveQuadraticSubMesh ) {
-      quadraticSubMeshAndViscousLayer = true;
-      _haveViscousLayers = !haveQuadraticSubMesh;
-      _comment += "Warning: Viscous layer is not possible with a quadratic mesh, it is ignored.";
-      error(COMPERR_WARNING, _comment);
-    }
-    // To remove as soon as quadratic mesh is allowed - END
-    
-    // needed to prevent the opencascade memory managmement from freeing things
-    vector<Handle(Geom2d_Curve)> curves;
-    vector<Handle(Geom_Surface)> surfaces;
-
-    fmap.Clear();
-    emap.Clear();
-    pmap.Clear();
-    FaceId2PythonSmp.clear();
-    EdgeId2PythonSmp.clear();
-    VertexId2PythonSmp.clear();
+  FacesWithSizeMap.Clear();
+  FaceId2SizeMap.clear();
+  FaceId2ClassAttractor.clear();
+  FaceIndex2ClassAttractor.clear();
+  EdgesWithSizeMap.Clear();
+  EdgeId2SizeMap.clear();
+  VerticesWithSizeMap.Clear();
+  VertexId2SizeMap.clear();
 
-    /****************************************************************************************
-                                          FACES
-    *****************************************************************************************/
-    int iface = 0;
-    string bad_end = "return";
-    int faceKey = -1;
-    TopTools_IndexedMapOfShape _map;
-    TopExp::MapShapes(aShape,TopAbs_VERTEX,_map);
-    int ienf = _map.Extent();
-
-    // Compute viscous layers at 2nd loop
-    {
-      TopTools_MapOfShape map;
-      if ( iLoop == 1 )
-        for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
-        {
-          const TopoDS_Face& f = TopoDS::Face(face_iter.Current());
-          if ( !map.Add( f )) continue;
-          iface = map.Extent();
-          viscousMeshPerFace[ iface ] = StdMeshers_ViscousLayers2D::Compute( aMesh, f );
-          if ( !viscousMeshPerFace[ iface ] )
-            return false; // error in StdMeshers_ViscousLayers2D::Compute()
-          if ( viscousMeshPerFace[ iface ]->NbProxySubMeshes() < 1 )
-            viscousMeshPerFace[ iface ].reset();
-        }
-      else
-      {
-        for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
-          map.Add( face_iter.Current());
-        viscousMeshPerFace.resize( map.Extent() + 1 );
-      }
-      iface = 0;
-    }
+  /* Now fill the CAD object with data from your CAD
+   * environement. This is the most complex part of a successfull
+   * integration.
+   */
 
-    assert(Py_IsInitialized());
-    PyGILState_STATE gstate;
-    gstate = PyGILState_Ensure();
+  // PreCAD
+  // If user requests it, send the CAD through Distene preprocessor : PreCAD
+  cad_t *cleanc = NULL; // preprocessed cad
+  precad_session_t *pcs = precad_session_new(ctx);
+  precad_data_set_cad(pcs, c);
+
+  cadsurf_session_t *css = cadsurf_session_new(ctx);
+
+  // an object that correctly deletes all cadsurf objects at destruction
+  BLSURF_Cleaner cleaner( ctx,css,c,dcad );
+
+  MESSAGE("BEGIN SetParameters");
+  bool use_precad = false;
+  SetParameters(
+                // #if BLSURF_VERSION_LONG >= "3.1.1"
+                //     c,
+                // #endif
+                _hypothesis, css, pcs, aMesh, &use_precad);
+  MESSAGE("END SetParameters");
+
+  haveQuadraticSubMesh = haveQuadraticSubMesh || (_hypothesis != NULL && _hypothesis->GetQuadraticMesh());
+  helper.SetIsQuadratic( haveQuadraticSubMesh );
+
+  // To remove as soon as quadratic mesh is allowed - BEGIN
+  // GDD: Viscous layer is not allowed with quadratic mesh
+  if (_haveViscousLayers && haveQuadraticSubMesh ) {
+    quadraticSubMeshAndViscousLayer = true;
+    _haveViscousLayers = !haveQuadraticSubMesh;
+    _comment += "Warning: Viscous layer is not possible with a quadratic mesh, it is ignored.";
+    error(COMPERR_WARNING, _comment);
+  }
+  // To remove as soon as quadratic mesh is allowed - END
 
-    string theSizeMapStr;
+  // needed to prevent the opencascade memory managmement from freeing things
+  vector<Handle(Geom2d_Curve)> curves;
+  vector<Handle(Geom_Surface)> surfaces;
 
-    for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
-    {
-      TopoDS_Face f = TopoDS::Face(face_iter.Current());
+  fmap.Clear();
+  emap.Clear();
+  pmap.Clear();
+  FaceId2PythonSmp.clear();
+  EdgeId2PythonSmp.clear();
+  VertexId2PythonSmp.clear();
 
-      // make INTERNAL face oriented FORWARD (issue 0020993)
-      if (f.Orientation() != TopAbs_FORWARD && f.Orientation() != TopAbs_REVERSED )
-        f.Orientation(TopAbs_FORWARD);
+  /****************************************************************************************
+                                          FACES
+  *****************************************************************************************/
+  int iface = 0;
+  string bad_end = "return";
+  int faceKey = -1;
+  TopTools_IndexedMapOfShape _map;
+  TopExp::MapShapes(aShape,TopAbs_VERTEX,_map);
+  int ienf = _map.Extent();
 
-      if (fmap.FindIndex(f) > 0)
-        continue;
-      iface = fmap.Add(f);
+  assert(Py_IsInitialized());
+  PyGILState_STATE gstate;
+  gstate = PyGILState_Ensure();
 
-      if ( iLoop == 1 && !viscousMeshPerFace[ iface ])
-        continue; // mesh already computed on f
+  string theSizeMapStr;
 
-      surfaces.push_back(BRep_Tool::Surface(f));
+  for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
+  {
+    TopoDS_Face f = TopoDS::Face(face_iter.Current());
 
-      /* create an object representing the face for cadsurf */
-      /* where face_id is an integer identifying the face.
-       * surf_function is the function that defines the surface
-       * (For this face, it will be called by cadsurf with your_face_object_ptr
-       * as last parameter.
-       */
-      cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
+    // make INTERNAL face oriented FORWARD (issue 0020993)
+    if (f.Orientation() != TopAbs_FORWARD && f.Orientation() != TopAbs_REVERSED )
+      f.Orientation(TopAbs_FORWARD);
 
-      /* by default a face has no tag (color).
-         The following call sets it to the same value as the face_id : */
-      cad_face_set_tag(fce, iface);
+    if (fmap.FindIndex(f) > 0)
+      continue;
+    iface = fmap.Add(f);
 
-      /* Set face orientation (optional if you want a well oriented output mesh)*/
-      if(f.Orientation() != TopAbs_FORWARD)
-        cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
-      else 
-        cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
+    surfaces.push_back(BRep_Tool::Surface(f));
 
-      if (HasSizeMapOnFace && !use_precad)
-      {
-        // -----------------
-        // Classic size map
-        // -----------------
-        faceKey = FacesWithSizeMap.FindIndex(f);
+    /* create an object representing the face for cadsurf */
+    /* where face_id is an integer identifying the face.
+     * surf_function is the function that defines the surface
+     * (For this face, it will be called by cadsurf with your_face_object_ptr
+     * as last parameter.
+     */
+    cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 
+    /* by default a face has no tag (color).
+       The following call sets it to the same value as the face_id : */
+    cad_face_set_tag(fce, iface);
 
-        if (FaceId2SizeMap.find(faceKey)!=FaceId2SizeMap.end()) {
-          MESSAGE("A size map is defined on face :"<<faceKey)
-            theSizeMapStr = FaceId2SizeMap[faceKey];
-          // check if function ends with "return"
-          if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
-            continue;
-          // Expr To Python function, verification is performed at validation in GUI
-          PyObject * obj = NULL;
-          obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
-          Py_DECREF(obj);
-          PyObject * func = NULL;
-          func = PyObject_GetAttrString(main_mod, "f");
-          FaceId2PythonSmp[iface]=func;
-          FaceId2SizeMap.erase(faceKey);
-        }
+    /* Set face orientation (optional if you want a well oriented output mesh)*/
+    if(f.Orientation() != TopAbs_FORWARD)
+      cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
+    else
+      cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 
-        // Specific size map = Attractor
-        std::map<int,std::vector<double> >::iterator attractor_iter = FaceId2AttractorCoords.begin();
-
-        for (; attractor_iter != FaceId2AttractorCoords.end(); ++attractor_iter) {
-          if (attractor_iter->first == faceKey) {
-            MESSAGE("Face indice: " << iface);
-            MESSAGE("Adding attractor");
-
-            double xyzCoords[3]  = {attractor_iter->second[2],
-                                    attractor_iter->second[3],
-                                    attractor_iter->second[4]};
-
-            MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
-            gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
-            BRepClass_FaceClassifier scl(f,P,1e-7);
-            // OCC 6.3sp6 : scl.Perform() is bugged. The function was rewritten
-            // BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
-            // OCC 6.5.2: scl.Perform() is not bugged anymore
-            scl.Perform(f, P, 1e-7);
-            TopAbs_State result = scl.State();
-            MESSAGE("Position of point on face: "<<result);
-            if ( result == TopAbs_OUT )
-              MESSAGE("Point is out of face: node is not created");
-            if ( result == TopAbs_UNKNOWN )
-              MESSAGE("Point position on face is unknown: node is not created");
-            if ( result == TopAbs_ON )
-              MESSAGE("Point is on border of face: node is not created");
-            if ( result == TopAbs_IN )
-            {
-              // Point is inside face and not on border
-              MESSAGE("Point is in face: node is created");
-              double uvCoords[2] = {attractor_iter->second[0],attractor_iter->second[1]};
-              ienf++;
-              MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
-              cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
-              cad_point_set_tag(point_p, ienf);
-            }
-            FaceId2AttractorCoords.erase(faceKey);
-          }
-        }
+    if (HasSizeMapOnFace && !use_precad)
+    {
+      // -----------------
+      // Classic size map
+      // -----------------
+      faceKey = FacesWithSizeMap.FindIndex(f);
 
-        // -----------------
-        // Class Attractors
-        // -----------------
-        std::map<int,BLSURFPlugin_Attractor* >::iterator clAttractor_iter = FaceId2ClassAttractor.find(faceKey);
-        if (clAttractor_iter != FaceId2ClassAttractor.end()){
+
+      if (FaceId2SizeMap.find(faceKey)!=FaceId2SizeMap.end()) {
+        MESSAGE("A size map is defined on face :"<<faceKey)
+          theSizeMapStr = FaceId2SizeMap[faceKey];
+        // check if function ends with "return"
+        if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
+          continue;
+        // Expr To Python function, verification is performed at validation in GUI
+        PyObject * obj = NULL;
+        obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
+        Py_DECREF(obj);
+        PyObject * func = NULL;
+        func = PyObject_GetAttrString(main_mod, "f");
+        FaceId2PythonSmp[iface]=func;
+        FaceId2SizeMap.erase(faceKey);
+      }
+
+      // Specific size map = Attractor
+      std::map<int,std::vector<double> >::iterator attractor_iter = FaceId2AttractorCoords.begin();
+
+      for (; attractor_iter != FaceId2AttractorCoords.end(); ++attractor_iter) {
+        if (attractor_iter->first == faceKey) {
           MESSAGE("Face indice: " << iface);
           MESSAGE("Adding attractor");
-          FaceIndex2ClassAttractor[iface]=clAttractor_iter->second;
-          FaceId2ClassAttractor.erase(clAttractor_iter);
-        }
-      } // if (HasSizeMapOnFace && !use_precad)
 
-      // ------------------
-      // Enforced Vertices
-      // ------------------
-      faceKey = FacesWithEnforcedVertices.FindIndex(f);
-      std::map<int,BLSURFPlugin_Hypothesis::TEnfVertexCoordsList >::const_iterator evmIt = FaceId2EnforcedVertexCoords.find(faceKey);
-      if (evmIt != FaceId2EnforcedVertexCoords.end()) {
-        MESSAGE("Some enforced vertices are defined");
-        BLSURFPlugin_Hypothesis::TEnfVertexCoordsList evl;
-        MESSAGE("Face indice: " << iface);
-        MESSAGE("Adding enforced vertices");
-        evl = evmIt->second;
-        MESSAGE("Number of vertices to add: "<< evl.size());
-        BLSURFPlugin_Hypothesis::TEnfVertexCoordsList::const_iterator evlIt = evl.begin();
-        for (; evlIt != evl.end(); ++evlIt) {
-          BLSURFPlugin_Hypothesis::TEnfVertexCoords xyzCoords;
-          xyzCoords.push_back(evlIt->at(2));
-          xyzCoords.push_back(evlIt->at(3));
-          xyzCoords.push_back(evlIt->at(4));
+          double xyzCoords[3]  = {attractor_iter->second[2],
+                                  attractor_iter->second[3],
+                                  attractor_iter->second[4]};
+
           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
           BRepClass_FaceClassifier scl(f,P,1e-7);
@@ -1704,694 +1756,700 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
           scl.Perform(f, P, 1e-7);
           TopAbs_State result = scl.State();
           MESSAGE("Position of point on face: "<<result);
-          if ( result == TopAbs_OUT ) {
+          if ( result == TopAbs_OUT )
             MESSAGE("Point is out of face: node is not created");
-            if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
-              EnfVertexCoords2ProjVertex.erase(xyzCoords);
-              EnfVertexCoords2EnfVertexList.erase(xyzCoords);
-            }
-          }
-          if ( result == TopAbs_UNKNOWN ) {
+          if ( result == TopAbs_UNKNOWN )
             MESSAGE("Point position on face is unknown: node is not created");
-            if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
-              EnfVertexCoords2ProjVertex.erase(xyzCoords);
-              EnfVertexCoords2EnfVertexList.erase(xyzCoords);
-            }
-          }
-          if ( result == TopAbs_ON ) {
+          if ( result == TopAbs_ON )
             MESSAGE("Point is on border of face: node is not created");
-            if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
-              EnfVertexCoords2ProjVertex.erase(xyzCoords);
-              EnfVertexCoords2EnfVertexList.erase(xyzCoords);
-            }
-          }
           if ( result == TopAbs_IN )
           {
             // Point is inside face and not on border
             MESSAGE("Point is in face: node is created");
-            double uvCoords[2]   = {evlIt->at(0),evlIt->at(1)};
+            double uvCoords[2] = {attractor_iter->second[0],attractor_iter->second[1]};
             ienf++;
             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
             cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
-            int tag = 0;
-            std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(xyzCoords);
-            if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end() &&
-                !enfCoordsIt->second.empty() )
-            {
-              // to merge nodes of an INTERNAL vertex belonging to several faces
-              TopoDS_Vertex     v = (*enfCoordsIt->second.begin())->vertex;
-              if ( v.IsNull() ) v = (*enfCoordsIt->second.rbegin())->vertex;
-              if ( !v.IsNull() ) {
-                tag = pmap.Add( v );
-                SMESH_subMesh* vSM = aMesh.GetSubMesh( v );
-                vSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );
-                mergeSubmeshes.insert( vSM->GetSubMeshDS() );
-                //if ( tag != pmap.Extent() )
-                  needMerge = true;
-              }
-            }
-            if ( tag == 0 ) tag = ienf;
-            cad_point_set_tag(point_p, tag);
+            cad_point_set_tag(point_p, ienf);
           }
+          FaceId2AttractorCoords.erase(faceKey);
         }
-        FaceId2EnforcedVertexCoords.erase(faceKey);
+      }
 
-      } // evmIt != FaceId2EnforcedVertexCoords.end()
-      //       else
-      //         std::cout << "No enforced vertex defined" << std::endl;
-      //     }
+      // -----------------
+      // Class Attractors
+      // -----------------
+      std::map<int,BLSURFPlugin_Attractor* >::iterator clAttractor_iter = FaceId2ClassAttractor.find(faceKey);
+      if (clAttractor_iter != FaceId2ClassAttractor.end()){
+        MESSAGE("Face indice: " << iface);
+        MESSAGE("Adding attractor");
+        FaceIndex2ClassAttractor[iface]=clAttractor_iter->second;
+        FaceId2ClassAttractor.erase(clAttractor_iter);
+      }
+    } // if (HasSizeMapOnFace && !use_precad)
 
-      /****************************************************************************************
-                                           EDGES
-                        now create the edges associated to this face
-      *****************************************************************************************/
-      int edgeKey = -1;
-      for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next())
-      {
-        TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
-        int ic = emap.FindIndex(e);
-        if (ic <= 0)
-          ic = emap.Add(e);
-
-        const SMESH_ProxyMesh::SubMesh* proxySubMesh = 0;
-        if ( viscousMeshPerFace[ iface ])
-          proxySubMesh = viscousMeshPerFace[ iface ]->GetProxySubMesh( e );
-
-        double tmin,tmax;
-        curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
-        if ( proxySubMesh )
-          curves.back() = new TProxyPCurve( proxySubMesh, curves.back() );
-
-        if (HasSizeMapOnEdge){
-          edgeKey = EdgesWithSizeMap.FindIndex(e);
-          if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end()) {
-            MESSAGE("Sizemap defined on edge with index " << edgeKey);
-            theSizeMapStr = EdgeId2SizeMap[edgeKey];
-            if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
-              continue;
-            // Expr To Python function, verification is performed at validation in GUI
-            PyObject * obj = NULL;
-            obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
-            Py_DECREF(obj);
-            PyObject * func = NULL;
-            func = PyObject_GetAttrString(main_mod, "f");
-            EdgeId2PythonSmp[ic]=func;
-            EdgeId2SizeMap.erase(edgeKey);
+      // ------------------
+      // Enforced Vertices
+      // ------------------
+    faceKey = FacesWithEnforcedVertices.FindIndex(f);
+    std::map<int,BLSURFPlugin_Hypothesis::TEnfVertexCoordsList >::const_iterator evmIt = FaceId2EnforcedVertexCoords.find(faceKey);
+    if (evmIt != FaceId2EnforcedVertexCoords.end()) {
+      MESSAGE("Some enforced vertices are defined");
+      BLSURFPlugin_Hypothesis::TEnfVertexCoordsList evl;
+      MESSAGE("Face indice: " << iface);
+      MESSAGE("Adding enforced vertices");
+      evl = evmIt->second;
+      MESSAGE("Number of vertices to add: "<< evl.size());
+      BLSURFPlugin_Hypothesis::TEnfVertexCoordsList::const_iterator evlIt = evl.begin();
+      for (; evlIt != evl.end(); ++evlIt) {
+        BLSURFPlugin_Hypothesis::TEnfVertexCoords xyzCoords;
+        xyzCoords.push_back(evlIt->at(2));
+        xyzCoords.push_back(evlIt->at(3));
+        xyzCoords.push_back(evlIt->at(4));
+        MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
+        gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
+        BRepClass_FaceClassifier scl(f,P,1e-7);
+        // OCC 6.3sp6 : scl.Perform() is bugged. The function was rewritten
+        // BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
+        // OCC 6.5.2: scl.Perform() is not bugged anymore
+        scl.Perform(f, P, 1e-7);
+        TopAbs_State result = scl.State();
+        MESSAGE("Position of point on face: "<<result);
+        if ( result == TopAbs_OUT ) {
+          MESSAGE("Point is out of face: node is not created");
+          if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
+            EnfVertexCoords2ProjVertex.erase(xyzCoords);
+            EnfVertexCoords2EnfVertexList.erase(xyzCoords);
           }
         }
-        /* data of nodes existing on the edge */
-        StdMeshers_FaceSidePtr nodeData;
-        SMESH_subMesh* sm = aMesh.GetSubMesh( e );
-        if ( !sm->IsEmpty() )
+        if ( result == TopAbs_UNKNOWN ) {
+          MESSAGE("Point position on face is unknown: node is not created");
+          if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
+            EnfVertexCoords2ProjVertex.erase(xyzCoords);
+            EnfVertexCoords2EnfVertexList.erase(xyzCoords);
+          }
+        }
+        if ( result == TopAbs_ON ) {
+          MESSAGE("Point is on border of face: node is not created");
+          if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
+            EnfVertexCoords2ProjVertex.erase(xyzCoords);
+            EnfVertexCoords2EnfVertexList.erase(xyzCoords);
+          }
+        }
+        if ( result == TopAbs_IN )
         {
-          SMESH_subMeshIteratorPtr subsmIt = sm->getDependsOnIterator( /*includeSelf=*/true,
-                                                                       /*complexFirst=*/false);
-          while ( subsmIt->more() ) 
-            edgeSubmeshes.insert( subsmIt->next()->GetSubMeshDS() );
-
-          if ( viscousMeshPerFace[ iface ])
-            if ( const SMESHDS_SubMesh* smDS = viscousMeshPerFace[ iface ]->GetProxySubMesh( e ))
-            { // protect proxySubmeshes from clearing
-              proxySubmeshes.insert( const_cast<SMESHDS_SubMesh*>( smDS ));
-              proxySubmeshes.insert( meshDS->MeshElements( e ));
-              proxySubmeshes.insert( meshDS->MeshElements( helper.IthVertex( 0, e )));
-              proxySubmeshes.insert( meshDS->MeshElements( helper.IthVertex( 1, e )));
-            }
-          nodeData.reset( new StdMeshers_FaceSide( f, e, &aMesh, /*isForwrd = */true,
-                                                   /*ignoreMedium=*/haveQuadraticSubMesh,
-                                                   viscousMeshPerFace[ iface ]));
-          if ( nodeData->MissVertexNode() )
-            return error(COMPERR_BAD_INPUT_MESH,"No node on vertex");
-
-          const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
-          if ( !nodeDataVec.empty() )
+          // Point is inside face and not on border
+          MESSAGE("Point is in face: node is created");
+          double uvCoords[2]   = {evlIt->at(0),evlIt->at(1)};
+          ienf++;
+          MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
+          cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
+          int tag = 0;
+          std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(xyzCoords);
+          if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end() &&
+              !enfCoordsIt->second.empty() )
           {
-            if ( Abs( nodeDataVec[0].param - tmin ) > Abs( nodeDataVec.back().param - tmin ))
-            {
-              nodeData->Reverse();
-              nodeData->GetUVPtStruct(); // nodeData recomputes nodeDataVec
+            // to merge nodes of an INTERNAL vertex belonging to several faces
+            TopoDS_Vertex     v = (*enfCoordsIt->second.begin())->vertex;
+            if ( v.IsNull() ) v = (*enfCoordsIt->second.rbegin())->vertex;
+            if ( !v.IsNull() ) {
+              tag = pmap.Add( v );
+              SMESH_subMesh* vSM = aMesh.GetSubMesh( v );
+              vSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );
+              mergeSubmeshes.insert( vSM->GetSubMeshDS() );
+              //if ( tag != pmap.Extent() )
+              needMerge = true;
             }
-            // tmin and tmax can change in case of viscous layer on an adjacent edge
-            tmin = nodeDataVec.front().param;
-            tmax = nodeDataVec.back().param;
-            //cout << "tmin = "<<tmin << ", tmax = "<< tmax << endl;
           }
-          else
+          if ( tag == 0 ) tag = ienf;
+          cad_point_set_tag(point_p, tag);
+        }
+      }
+      FaceId2EnforcedVertexCoords.erase(faceKey);
+
+    }
+
+    /****************************************************************************************
+                                           EDGES
+                        now create the edges associated to this face
+    *****************************************************************************************/
+    int edgeKey = -1;
+    for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next())
+    {
+      TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
+      int ic = emap.FindIndex(e);
+      if (ic <= 0)
+        ic = emap.Add(e);
+
+      double tmin,tmax;
+      curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
+
+      if (HasSizeMapOnEdge){
+        edgeKey = EdgesWithSizeMap.FindIndex(e);
+        if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end()) {
+          MESSAGE("Sizemap defined on edge with index " << edgeKey);
+          theSizeMapStr = EdgeId2SizeMap[edgeKey];
+          if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
+            continue;
+          // Expr To Python function, verification is performed at validation in GUI
+          PyObject * obj = NULL;
+          obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
+          Py_DECREF(obj);
+          PyObject * func = NULL;
+          func = PyObject_GetAttrString(main_mod, "f");
+          EdgeId2PythonSmp[ic]=func;
+          EdgeId2SizeMap.erase(edgeKey);
+        }
+      }
+      /* data of nodes existing on the edge */
+      StdMeshers_FaceSidePtr nodeData;
+      SMESH_subMesh* sm = aMesh.GetSubMesh( e );
+      if ( !sm->IsEmpty() )
+      {
+        SMESH_subMeshIteratorPtr subsmIt = sm->getDependsOnIterator( /*includeSelf=*/true,
+                                                                     /*complexFirst=*/false);
+        while ( subsmIt->more() )
+          edgeSubmeshes.insert( subsmIt->next()->GetSubMeshDS() );
+
+        nodeData.reset( new StdMeshers_FaceSide( f, e, &aMesh, /*isForwrd = */true,
+                                                 /*ignoreMedium=*/haveQuadraticSubMesh));
+        if ( nodeData->MissVertexNode() )
+          return error(COMPERR_BAD_INPUT_MESH,"No node on vertex");
+
+        const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
+        if ( !nodeDataVec.empty() )
+        {
+          if ( Abs( nodeDataVec[0].param - tmin ) > Abs( nodeDataVec.back().param - tmin ))
           {
-            cout << "---------------- Invalid nodeData" << endl;
-            nodeData.reset();
+            nodeData->Reverse();
+            nodeData->GetUVPtStruct(); // nodeData recomputes nodeDataVec
           }
+          // tmin and tmax can change in case of viscous layer on an adjacent edge
+          tmin = nodeDataVec.front().param;
+          tmax = nodeDataVec.back().param;
+        }
+        else
+        {
+          cout << "---------------- Invalid nodeData" << endl;
+          nodeData.reset();
         }
+      }
 
-        /* attach the edge to the current cadsurf face */
-        cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
+      /* attach the edge to the current cadsurf face */
+      cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
 
-        /* by default an edge has no tag (color).
-           The following call sets it to the same value as the edge_id : */
-        cad_edge_set_tag(edg, ic);
+      /* by default an edge has no tag (color).
+         The following call sets it to the same value as the edge_id : */
+      cad_edge_set_tag(edg, ic);
 
-        /* by default, an edge does not necessalry appear in the resulting mesh,
-           unless the following property is set :
-        */
-        cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
+      /* by default, an edge does not necessalry appear in the resulting mesh,
+         unless the following property is set :
+      */
+      cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
 
-        /* by default an edge is a boundary edge */
-        if (e.Orientation() == TopAbs_INTERNAL)
-          cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
+      /* by default an edge is a boundary edge */
+      if (e.Orientation() == TopAbs_INTERNAL)
+        cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
 
-        // pass existing nodes of sub-meshes to BLSURF
-        if ( nodeData )
-        {
-          const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
-          const int                      nbNodes     = nodeDataVec.size();
+      // pass existing nodes of sub-meshes to BLSURF
+      if ( nodeData )
+      {
+        const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
+        const int                      nbNodes     = nodeDataVec.size();
 
-          dcad_edge_discretization_t *dedge;
-          dcad_get_edge_discretization(dcad, edg, &dedge);
-          dcad_edge_discretization_set_vertex_count( dedge, nbNodes );
+        dcad_edge_discretization_t *dedge;
+        dcad_get_edge_discretization(dcad, edg, &dedge);
+        dcad_edge_discretization_set_vertex_count( dedge, nbNodes );
 
-          //cout << endl << " EDGE " << ic << endl;
-          for ( int iN = 0; iN < nbNodes; ++iN )
-          {
-            const UVPtStruct& nData = nodeDataVec[ iN ];
-            double t                = nData.param;
-            real uv[2]              = { nData.u, nData.v };
-            SMESH_TNodeXYZ nXYZ( nData.node );
-            //cout << "\tt = " << t << " uv = ( " << uv[0] << ","<< uv[1] << " ) ID " << nData.node->GetID() << endl;
-            dcad_edge_discretization_set_vertex_coordinates( dedge, iN+1, t, uv, nXYZ._xyz );
-          }
-          dcad_edge_discretization_set_property(dedge, DISTENE_DCAD_PROPERTY_REQUIRED);
+        // cout << endl << " EDGE " << ic << endl;
+        // cout << "tmin = "<<tmin << ", tmax = "<< tmax << endl;
+        for ( int iN = 0; iN < nbNodes; ++iN )
+        {
+          const UVPtStruct& nData = nodeDataVec[ iN ];
+          double t                = nData.param;
+          real uv[2]              = { nData.u, nData.v };
+          SMESH_TNodeXYZ nXYZ( nData.node );
+          // cout << "\tt = " << t
+          //      << "\t uv = ( " << uv[0] << ","<< uv[1] << " ) "
+          //      << "\t u = " << nData.param
+          //      << "\t ID = " << nData.node->GetID() << endl;
+          dcad_edge_discretization_set_vertex_coordinates( dedge, iN+1, t, uv, nXYZ._xyz );
         }
+        dcad_edge_discretization_set_property(dedge, DISTENE_DCAD_PROPERTY_REQUIRED);
+      }
 
-        /****************************************************************************************
+      /****************************************************************************************
                                       VERTICES
-        *****************************************************************************************/
-
-        int npts = 0;
-        int ip1, ip2, *ip;
-        gp_Pnt2d e0 = curves.back()->Value(tmin);
-        gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
-        Standard_Real d1=0,d2=0;
-
-        int vertexKey = -1;
-        for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
-          TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
-          ++npts;
-          if (npts == 1){
-            ip = &ip1;
-            d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
-          } else {
-            ip = &ip2;
-            d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
-          }
-          *ip = pmap.FindIndex(v);
-          if(*ip <= 0)
-            *ip = pmap.Add(v);
-
-          if (HasSizeMapOnVertex){
-            vertexKey = VerticesWithSizeMap.FindIndex(v);
-            if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
-              theSizeMapStr = VertexId2SizeMap[vertexKey];
-              //MESSAGE("VertexId2SizeMap[faceKey]: " << VertexId2SizeMap[vertexKey]);
-              if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
-                continue;
-              // Expr To Python function, verification is performed at validation in GUI
-              PyObject * obj = NULL;
-              obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
-              Py_DECREF(obj);
-              PyObject * func = NULL;
-              func = PyObject_GetAttrString(main_mod, "f");
-              VertexId2PythonSmp[*ip]=func;
-              VertexId2SizeMap.erase(vertexKey);   // do not erase if using a vector
-            }
-          }
-        }
-        if (npts != 2) {
-          // should not happen
-          MESSAGE("An edge does not have 2 extremities.");
+      *****************************************************************************************/
+
+      int npts = 0;
+      int ip1, ip2, *ip;
+      gp_Pnt2d e0 = curves.back()->Value(tmin);
+      gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
+      Standard_Real d1=0,d2=0;
+
+      int vertexKey = -1;
+      for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
+        TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
+        ++npts;
+        if (npts == 1){
+          ip = &ip1;
+          d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
         } else {
-          if (d1 < d2) {
-            // This defines the curves extremity connectivity
-            cad_edge_set_extremities(edg, ip1, ip2);
-            /* set the tag (color) to the same value as the extremity id : */
-            cad_edge_set_extremities_tag(edg, ip1, ip2);
-          }
-          else {
-            cad_edge_set_extremities(edg, ip2, ip1);
-            cad_edge_set_extremities_tag(edg, ip2, ip1);
-          }
+          ip = &ip2;
+          d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
         }
-      } // for edge
-    } //for face
-
-    // Clear mesh from already meshed edges if possible else
-    // remember that merge is needed
-    set< SMESHDS_SubMesh* >::iterator smIt = edgeSubmeshes.begin();
-    for ( ; smIt != edgeSubmeshes.end(); ++smIt ) // loop on already meshed EDGEs
-    {
-      SMESHDS_SubMesh* smDS = *smIt;
-      if ( !smDS ) continue;
-      SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
-      if ( nIt->more() )
-      {
-        const SMDS_MeshNode* n = nIt->next();
-        if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
-        {
-          needMerge = true;
-          // add existing medium nodes to helper
-          if ( aMesh.NbEdges( ORDER_QUADRATIC ) > 0 )
-          {
-            SMDS_ElemIteratorPtr edgeIt = smDS->GetElements();
-            while ( edgeIt->more() )
-              helper.AddTLinks( static_cast<const SMDS_MeshEdge*>(edgeIt->next()));
+        *ip = pmap.FindIndex(v);
+        if(*ip <= 0) {
+          *ip = pmap.Add(v);
+          SMESH_subMesh* sm = aMesh.GetSubMesh(v);
+          if ( sm->IsMeshComputed() )
+            edgeSubmeshes.insert( sm->GetSubMeshDS() );
+        }
+        if (HasSizeMapOnVertex){
+          vertexKey = VerticesWithSizeMap.FindIndex(v);
+          if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
+            theSizeMapStr = VertexId2SizeMap[vertexKey];
+            //MESSAGE("VertexId2SizeMap[faceKey]: " << VertexId2SizeMap[vertexKey]);
+            if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
+              continue;
+            // Expr To Python function, verification is performed at validation in GUI
+            PyObject * obj = NULL;
+            obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
+            Py_DECREF(obj);
+            PyObject * func = NULL;
+            func = PyObject_GetAttrString(main_mod, "f");
+            VertexId2PythonSmp[*ip]=func;
+            VertexId2SizeMap.erase(vertexKey);   // do not erase if using a vector
           }
-          continue;
         }
       }
-      if ( proxySubmeshes.count( smDS ))
-      {
-        needMerge = true;
-      }
-      else
-      {
-        SMDS_ElemIteratorPtr eIt = smDS->GetElements();
-        while ( eIt->more() ) meshDS->RemoveFreeElement( eIt->next(), smDS );
-        SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
-        while ( nIt->more() ) meshDS->RemoveFreeNode( nIt->next(), smDS );
+      if (npts != 2) {
+        // should not happen
+        MESSAGE("An edge does not have 2 extremities.");
+      } else {
+        if (d1 < d2) {
+          // This defines the curves extremity connectivity
+          cad_edge_set_extremities(edg, ip1, ip2);
+          /* set the tag (color) to the same value as the extremity id : */
+          cad_edge_set_extremities_tag(edg, ip1, ip2);
+        }
+        else {
+          cad_edge_set_extremities(edg, ip2, ip1);
+          cad_edge_set_extremities_tag(edg, ip2, ip1);
+        }
       }
-      if ( proxySubmeshes.count( smDS ))
+    } // for edge
+  } //for face
+
+  // Clear mesh from already meshed edges if possible else
+  // remember that merge is needed
+  TSubMeshSet::iterator smIt = edgeSubmeshes.begin();
+  for ( ; smIt != edgeSubmeshes.end(); ++smIt ) // loop on already meshed EDGEs
+  {
+    SMESHDS_SubMesh* smDS = *smIt;
+    if ( !smDS ) continue;
+    SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
+    if ( nIt->more() )
+    {
+      const SMDS_MeshNode* n = nIt->next();
+      if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
       {
         needMerge = true;
-      }
-      else
-      {
-        SMDS_ElemIteratorPtr eIt = smDS->GetElements();
-        while ( eIt->more() ) meshDS->RemoveFreeElement( eIt->next(), smDS );
-        SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
-        while ( nIt->more() ) meshDS->RemoveFreeNode( nIt->next(), smDS );
-      }
-    }
-
-    if (use_precad) {
-      /* Now launch the PreCAD process */
-      status = precad_process(pcs);
-      if(status != STATUS_OK){
-        cout << "PreCAD processing failed with error code " << status << "\n";
-      }
-      else {
-        // retrieve the pre-processed CAD object
-        cleanc = precad_new_cad(pcs);
-        if(!cleanc){
-          cout << "Unable to retrieve PreCAD result \n";
+        // add existing medium nodes to helper
+        if ( aMesh.NbEdges( ORDER_QUADRATIC ) > 0 )
+        {
+          SMDS_ElemIteratorPtr edgeIt = smDS->GetElements();
+          while ( edgeIt->more() )
+            helper.AddTLinks( static_cast<const SMDS_MeshEdge*>(edgeIt->next()));
         }
-        cout << "PreCAD processing successfull \n";
-
-        // #if BLSURF_VERSION_LONG >= "3.1.1"
-        //       /* We can now get the updated sizemaps (if any) */
-        // //       if(geo_sizemap_e)
-        // //         clean_geo_sizemap_e = precad_new_sizemap(pcs, geo_sizemap_e);
-        // // 
-        // //       if(geo_sizemap_f)
-        // //         clean_geo_sizemap_f = precad_new_sizemap(pcs, geo_sizemap_f);
-        //
-        //       if(iso_sizemap_p)
-        //         clean_iso_sizemap_p = precad_new_sizemap(pcs, iso_sizemap_p);
-        //
-        //       if(iso_sizemap_e)
-        //         clean_iso_sizemap_e = precad_new_sizemap(pcs, iso_sizemap_e);
-        //
-        //       if(iso_sizemap_f)
-        //         clean_iso_sizemap_f = precad_new_sizemap(pcs, iso_sizemap_f);
-        // #endif
+        continue;
       }
-      // Now we can delete the PreCAD session
-      precad_session_delete(pcs);
     }
+    {
+      SMDS_ElemIteratorPtr eIt = smDS->GetElements();
+      while ( eIt->more() ) meshDS->RemoveFreeElement( eIt->next(), smDS );
+      SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
+      while ( nIt->more() ) meshDS->RemoveFreeNode( nIt->next(), smDS );
+    }
+  }
 
-    cadsurf_data_set_dcad(css, dcad);
-    if (cleanc) {
-      // Give the pre-processed CAD object to the current BLSurf session
-      cadsurf_data_set_cad(css, cleanc);
+
+  if (use_precad) {
+    /* Now launch the PreCAD process */
+    status = precad_process(pcs);
+    if(status != STATUS_OK){
+      cout << "PreCAD processing failed with error code " << status << "\n";
     }
     else {
-      // Use the original one
-      cadsurf_data_set_cad(css, c);
+      // retrieve the pre-processed CAD object
+      cleanc = precad_new_cad(pcs);
+      if(!cleanc){
+        cout << "Unable to retrieve PreCAD result \n";
+      }
+      cout << "PreCAD processing successfull \n";
+
+      // #if BLSURF_VERSION_LONG >= "3.1.1"
+      //       /* We can now get the updated sizemaps (if any) */
+      // //       if(geo_sizemap_e)
+      // //         clean_geo_sizemap_e = precad_new_sizemap(pcs, geo_sizemap_e);
+      // // 
+      // //       if(geo_sizemap_f)
+      // //         clean_geo_sizemap_f = precad_new_sizemap(pcs, geo_sizemap_f);
+      //
+      //       if(iso_sizemap_p)
+      //         clean_iso_sizemap_p = precad_new_sizemap(pcs, iso_sizemap_p);
+      //
+      //       if(iso_sizemap_e)
+      //         clean_iso_sizemap_e = precad_new_sizemap(pcs, iso_sizemap_e);
+      //
+      //       if(iso_sizemap_f)
+      //         clean_iso_sizemap_f = precad_new_sizemap(pcs, iso_sizemap_f);
+      // #endif
     }
+    // Now we can delete the PreCAD session
+    precad_session_delete(pcs);
+  }
 
-    std::cout << std::endl;
-    std::cout << "Beginning of Surface Mesh generation" << std::endl;
-    std::cout << std::endl;
+  cadsurf_data_set_dcad(css, dcad);
+  if (cleanc) {
+    // Give the pre-processed CAD object to the current BLSurf session
+    cadsurf_data_set_cad(css, cleanc);
+  }
+  else {
+    // Use the original one
+    cadsurf_data_set_cad(css, c);
+  }
 
-    try {
-      OCC_CATCH_SIGNALS;
+  std::cout << std::endl;
+  std::cout << "Beginning of Surface Mesh generation" << std::endl;
+  std::cout << std::endl;
 
-      status = cadsurf_compute_mesh(css);
+  try {
+    OCC_CATCH_SIGNALS;
 
+    status = cadsurf_compute_mesh(css);
+
+  }
+  catch ( std::exception& exc ) {
+    _comment += exc.what();
+  }
+  catch (Standard_Failure& ex) {
+    _comment += ex.DynamicType()->Name();
+    if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
+      _comment += ": ";
+      _comment += ex.GetMessageString();
     }
-    catch ( std::exception& exc ) {
-      _comment += exc.what();
-    }
-    catch (Standard_Failure& ex) {
-      _comment += ex.DynamicType()->Name();
-      if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
-        _comment += ": ";
-        _comment += ex.GetMessageString();
-      }
-    }
-    catch (...) {
-      if ( _comment.empty() )
-        _comment = "Exception in cadsurf_compute_mesh()";
-    }
-    if ( status != STATUS_OK) {
-      // There was an error while meshing
-      error(_comment);
-    }
+  }
+  catch (...) {
+    if ( _comment.empty() )
+      _comment = "Exception in cadsurf_compute_mesh()";
+  }
+  if ( status != STATUS_OK) {
+    // There was an error while meshing
+    error(_comment);
+  }
 
-    PyGILState_Release(gstate);
+  PyGILState_Release(gstate);
 
-    std::cout << std::endl;
-    std::cout << "End of Surface Mesh generation" << std::endl;
-    std::cout << std::endl;
+  std::cout << std::endl;
+  std::cout << "End of Surface Mesh generation" << std::endl;
+  std::cout << std::endl;
 
-    mesh_t *msh = NULL;
-    cadsurf_data_get_mesh(css, &msh);
-    if(!msh){
-      /* release the mesh object */
-      cadsurf_data_regain_mesh(css, msh);
-      return error(_comment);
-    }
+  mesh_t *msh = NULL;
+  cadsurf_data_get_mesh(css, &msh);
+  if(!msh){
+    /* release the mesh object */
+    cadsurf_data_regain_mesh(css, msh);
+    return error(_comment);
+  }
 
-    std::string GMFFileName = BLSURFPlugin_Hypothesis::GetDefaultGMFFile();
-    if (_hypothesis)
-      GMFFileName = _hypothesis->GetGMFFile();
-    if (GMFFileName != "") {
-      //     bool GMFFileMode = _hypothesis->GetGMFFileMode();
-      bool asciiFound = (GMFFileName.find(".mesh",GMFFileName.length()-5) != std::string::npos);
-      bool binaryFound = (GMFFileName.find(".meshb",GMFFileName.length()-6) != std::string::npos);
-      if (!asciiFound && !binaryFound)
-        GMFFileName.append(".mesh");
-      mesh_write_mesh(msh, GMFFileName.c_str());
-    }
+  std::string GMFFileName = BLSURFPlugin_Hypothesis::GetDefaultGMFFile();
+  if (_hypothesis)
+    GMFFileName = _hypothesis->GetGMFFile();
+  if (GMFFileName != "") {
+    //     bool GMFFileMode = _hypothesis->GetGMFFileMode();
+    bool asciiFound = (GMFFileName.find(".mesh",GMFFileName.length()-5) != std::string::npos);
+    bool binaryFound = (GMFFileName.find(".meshb",GMFFileName.length()-6) != std::string::npos);
+    if (!asciiFound && !binaryFound)
+      GMFFileName.append(".mesh");
+    mesh_write_mesh(msh, GMFFileName.c_str());
+  }
 
-    /* retrieve mesh data (see meshgems/mesh.h) */
-    integer nv, ne, nt, nq, vtx[4], tag;
-    integer *evedg, *evtri, *evquad, type;
-    real xyz[3];
-
-    mesh_get_vertex_count(msh, &nv);
-    mesh_get_edge_count(msh, &ne);
-    mesh_get_triangle_count(msh, &nt);
-    mesh_get_quadrangle_count(msh, &nq);
-
-    evedg  = (integer *)mesh_calloc_generic_buffer(msh);
-    evtri  = (integer *)mesh_calloc_generic_buffer(msh);
-    evquad = (integer *)mesh_calloc_generic_buffer(msh);
-
-    SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
-    bool* tags = new bool[nv+1];
-
-    /* enumerated vertices */
-    for(int iv=1;iv<=nv;iv++) {
-      mesh_get_vertex_coordinates(msh, iv, xyz);
-      mesh_get_vertex_tag(msh, iv, &tag);
-      // Issue 0020656. Use vertex coordinates
-      if ( tag > 0 && tag <= pmap.Extent() ) {
-        TopoDS_Vertex v = TopoDS::Vertex(pmap(tag));
-        double tol = BRep_Tool::Tolerance( v );
-        gp_Pnt p = BRep_Tool::Pnt( v );
-        if ( p.IsEqual( gp_Pnt( xyz[0], xyz[1], xyz[2]), 2*tol))
-          xyz[0] = p.X(), xyz[1] = p.Y(), xyz[2] = p.Z();
-        else
-          tag = 0; // enforced or attracted vertex
-      }
-      nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
-
-      // Create group of enforced vertices if requested
-      BLSURFPlugin_Hypothesis::TEnfVertexCoords projVertex;
-      projVertex.clear();
-      projVertex.push_back((double)xyz[0]);
-      projVertex.push_back((double)xyz[1]);
-      projVertex.push_back((double)xyz[2]);
-      std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(projVertex);
-      if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end()) {
-        MESSAGE("Found enforced vertex @ " << xyz[0] << ", " << xyz[1] << ", " << xyz[2]);
-        BLSURFPlugin_Hypothesis::TEnfVertexList::const_iterator enfListIt = enfCoordsIt->second.begin();
-        BLSURFPlugin_Hypothesis::TEnfVertex *currentEnfVertex;
-        for (; enfListIt != enfCoordsIt->second.end(); ++enfListIt) {
-          currentEnfVertex = (*enfListIt);
-          if (currentEnfVertex->grpName != "") {
-            bool groupDone = false;
-            SMESH_Mesh::GroupIteratorPtr grIt = aMesh.GetGroups();
-            MESSAGE("currentEnfVertex->grpName: " << currentEnfVertex->grpName);
-            MESSAGE("Parsing the groups of the mesh");
-            while (grIt->more()) {
-              SMESH_Group * group = grIt->next();
-              if ( !group ) continue;
-              MESSAGE("Group: " << group->GetName());
-              SMESHDS_GroupBase* groupDS = group->GetGroupDS();
-              if ( !groupDS ) continue;
-              MESSAGE("group->SMDSGroup().GetType(): " << (groupDS->GetType()));
-              MESSAGE("group->SMDSGroup().GetType()==SMDSAbs_Node: " << (groupDS->GetType()==SMDSAbs_Node));
-              MESSAGE("currentEnfVertex->grpName.compare(group->GetStoreName())==0: " << (currentEnfVertex->grpName.compare(group->GetName())==0));
-              if ( groupDS->GetType()==SMDSAbs_Node && currentEnfVertex->grpName.compare(group->GetName())==0) {
-                SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
-                aGroupDS->SMDSGroup().Add(nodes[iv]);
-                MESSAGE("Node ID: " << nodes[iv]->GetID());
-                // How can I inform the hypothesis ?
-                //                 _hypothesis->AddEnfVertexNodeID(currentEnfVertex->grpName,nodes[iv]->GetID());
-                groupDone = true;
-                MESSAGE("Successfully added enforced vertex to existing group " << currentEnfVertex->grpName);
-                break;
-              }
-            }
-            if (!groupDone)
-            {
-              int groupId;
-              SMESH_Group* aGroup = aMesh.AddGroup(SMDSAbs_Node, currentEnfVertex->grpName.c_str(), groupId);
-              aGroup->SetName( currentEnfVertex->grpName.c_str() );
-              SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
+  /* retrieve mesh data (see meshgems/mesh.h) */
+  integer nv, ne, nt, nq, vtx[4], tag;
+  integer *evedg, *evtri, *evquad, type;
+  real xyz[3];
+
+  mesh_get_vertex_count(msh, &nv);
+  mesh_get_edge_count(msh, &ne);
+  mesh_get_triangle_count(msh, &nt);
+  mesh_get_quadrangle_count(msh, &nq);
+
+  evedg  = (integer *)mesh_calloc_generic_buffer(msh);
+  evtri  = (integer *)mesh_calloc_generic_buffer(msh);
+  evquad = (integer *)mesh_calloc_generic_buffer(msh);
+
+  SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
+  bool* tags = new bool[nv+1];
+
+  /* enumerated vertices */
+  for(int iv=1;iv<=nv;iv++) {
+    mesh_get_vertex_coordinates(msh, iv, xyz);
+    mesh_get_vertex_tag(msh, iv, &tag);
+    // Issue 0020656. Use vertex coordinates
+    if ( tag > 0 && tag <= pmap.Extent() ) {
+      TopoDS_Vertex v = TopoDS::Vertex(pmap(tag));
+      double tol = BRep_Tool::Tolerance( v );
+      gp_Pnt p = BRep_Tool::Pnt( v );
+      if ( p.IsEqual( gp_Pnt( xyz[0], xyz[1], xyz[2]), 2*tol))
+        xyz[0] = p.X(), xyz[1] = p.Y(), xyz[2] = p.Z();
+      else
+        tag = 0; // enforced or attracted vertex
+    }
+    nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
+
+    // Create group of enforced vertices if requested
+    BLSURFPlugin_Hypothesis::TEnfVertexCoords projVertex;
+    projVertex.clear();
+    projVertex.push_back((double)xyz[0]);
+    projVertex.push_back((double)xyz[1]);
+    projVertex.push_back((double)xyz[2]);
+    std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(projVertex);
+    if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end()) {
+      MESSAGE("Found enforced vertex @ " << xyz[0] << ", " << xyz[1] << ", " << xyz[2]);
+      BLSURFPlugin_Hypothesis::TEnfVertexList::const_iterator enfListIt = enfCoordsIt->second.begin();
+      BLSURFPlugin_Hypothesis::TEnfVertex *currentEnfVertex;
+      for (; enfListIt != enfCoordsIt->second.end(); ++enfListIt) {
+        currentEnfVertex = (*enfListIt);
+        if (currentEnfVertex->grpName != "") {
+          bool groupDone = false;
+          SMESH_Mesh::GroupIteratorPtr grIt = aMesh.GetGroups();
+          MESSAGE("currentEnfVertex->grpName: " << currentEnfVertex->grpName);
+          MESSAGE("Parsing the groups of the mesh");
+          while (grIt->more()) {
+            SMESH_Group * group = grIt->next();
+            if ( !group ) continue;
+            MESSAGE("Group: " << group->GetName());
+            SMESHDS_GroupBase* groupDS = group->GetGroupDS();
+            if ( !groupDS ) continue;
+            MESSAGE("group->SMDSGroup().GetType(): " << (groupDS->GetType()));
+            MESSAGE("group->SMDSGroup().GetType()==SMDSAbs_Node: " << (groupDS->GetType()==SMDSAbs_Node));
+            MESSAGE("currentEnfVertex->grpName.compare(group->GetStoreName())==0: " << (currentEnfVertex->grpName.compare(group->GetName())==0));
+            if ( groupDS->GetType()==SMDSAbs_Node && currentEnfVertex->grpName.compare(group->GetName())==0) {
+              SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
               aGroupDS->SMDSGroup().Add(nodes[iv]);
-              MESSAGE("Successfully created enforced vertex group " << currentEnfVertex->grpName);
+              MESSAGE("Node ID: " << nodes[iv]->GetID());
+              // How can I inform the hypothesis ?
+              //                 _hypothesis->AddEnfVertexNodeID(currentEnfVertex->grpName,nodes[iv]->GetID());
               groupDone = true;
+              MESSAGE("Successfully added enforced vertex to existing group " << currentEnfVertex->grpName);
+              break;
             }
-            if (!groupDone)
-              throw SALOME_Exception(LOCALIZED("An enforced vertex node was not added to a group"));
           }
-          else
-            MESSAGE("Group name is empty: '"<<currentEnfVertex->grpName<<"' => group is not created");
+          if (!groupDone)
+          {
+            int groupId;
+            SMESH_Group* aGroup = aMesh.AddGroup(SMDSAbs_Node, currentEnfVertex->grpName.c_str(), groupId);
+            aGroup->SetName( currentEnfVertex->grpName.c_str() );
+            SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
+            aGroupDS->SMDSGroup().Add(nodes[iv]);
+            MESSAGE("Successfully created enforced vertex group " << currentEnfVertex->grpName);
+            groupDone = true;
+          }
+          if (!groupDone)
+            throw SALOME_Exception(LOCALIZED("An enforced vertex node was not added to a group"));
         }
+        else
+          MESSAGE("Group name is empty: '"<<currentEnfVertex->grpName<<"' => group is not created");
       }
+    }
 
-      // internal points are tagged to zero
-      if(tag > 0 && tag <= pmap.Extent() ){
-        meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
-        tags[iv] = false;
-      } else {
-        tags[iv] = true;
-      }
+    // internal points are tagged to zero
+    if(tag > 0 && tag <= pmap.Extent() ){
+      meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
+      tags[iv] = false;
+    } else {
+      tags[iv] = true;
     }
+  }
 
-    /* enumerate edges */
-    for(int it=1;it<=ne;it++) {
-      SMDS_MeshEdge* edg;
-      mesh_get_edge_vertices(msh, it, vtx);
-      mesh_get_edge_extra_vertices(msh, it, &type, evedg);
-      mesh_get_edge_tag(msh, it, &tag);
-      if (tags[vtx[0]]) {
-        Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
-        tags[vtx[0]] = false;
-      };
-      if (tags[vtx[1]]) {
-        Set_NodeOnEdge(meshDS, nodes[vtx[1]], emap(tag));
-        tags[vtx[1]] = false;
-      };
-      if (type == MESHGEMS_MESH_ELEMENT_TYPE_EDGE3) {
-        // QUADRATIC EDGE
-        if (tags[evedg[0]]) {
-          Set_NodeOnEdge(meshDS, nodes[evedg[0]], emap(tag));
-          tags[evedg[0]] = false;
-        }
-        edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]], nodes[evedg[0]]);
-      }
-      else {
-        edg = helper.AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
+  /* enumerate edges */
+  for(int it=1;it<=ne;it++) {
+    SMDS_MeshEdge* edg;
+    mesh_get_edge_vertices(msh, it, vtx);
+    mesh_get_edge_extra_vertices(msh, it, &type, evedg);
+    mesh_get_edge_tag(msh, it, &tag);
+    if (tags[vtx[0]]) {
+      Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
+      tags[vtx[0]] = false;
+    };
+    if (tags[vtx[1]]) {
+      Set_NodeOnEdge(meshDS, nodes[vtx[1]], emap(tag));
+      tags[vtx[1]] = false;
+    };
+    if (type == MESHGEMS_MESH_ELEMENT_TYPE_EDGE3) {
+      // QUADRATIC EDGE
+      if (tags[evedg[0]]) {
+        Set_NodeOnEdge(meshDS, nodes[evedg[0]], emap(tag));
+        tags[evedg[0]] = false;
       }
-      meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
+      edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]], nodes[evedg[0]]);
     }
+    else {
+      edg = helper.AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
+    }
+    meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
+  }
 
-    /* enumerate triangles */
-    for(int it=1;it<=nt;it++) {
-      SMDS_MeshFace* tri;
-      mesh_get_triangle_vertices(msh, it, vtx);
-      mesh_get_triangle_extra_vertices(msh, it, &type, evtri);
-      mesh_get_triangle_tag(msh, it, &tag);
-      if (tags[vtx[0]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
-        tags[vtx[0]] = false;
-      };
-      if (tags[vtx[1]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
-        tags[vtx[1]] = false;
-      };
-      if (tags[vtx[2]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
-        tags[vtx[2]] = false;
-      };
-      if (type == MESHGEMS_MESH_ELEMENT_TYPE_TRIA6) {
-        // QUADRATIC TRIANGLE
-        if (tags[evtri[0]]) {
-          meshDS->SetNodeOnFace(nodes[evtri[0]], TopoDS::Face(fmap(tag)));
-          tags[evtri[0]] = false;
-        }
-        if (tags[evtri[1]]) {
-          meshDS->SetNodeOnFace(nodes[evtri[1]], TopoDS::Face(fmap(tag)));
-          tags[evtri[1]] = false;
-        }
-        if (tags[evtri[2]]) {
-          meshDS->SetNodeOnFace(nodes[evtri[2]], TopoDS::Face(fmap(tag)));
-          tags[evtri[2]] = false;
-        }
-        tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]],
-                              nodes[evtri[0]], nodes[evtri[1]], nodes[evtri[2]]);
+  /* enumerate triangles */
+  for(int it=1;it<=nt;it++) {
+    SMDS_MeshFace* tri;
+    mesh_get_triangle_vertices(msh, it, vtx);
+    mesh_get_triangle_extra_vertices(msh, it, &type, evtri);
+    mesh_get_triangle_tag(msh, it, &tag);
+    if (tags[vtx[0]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
+      tags[vtx[0]] = false;
+    };
+    if (tags[vtx[1]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
+      tags[vtx[1]] = false;
+    };
+    if (tags[vtx[2]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
+      tags[vtx[2]] = false;
+    };
+    if (type == MESHGEMS_MESH_ELEMENT_TYPE_TRIA6) {
+      // QUADRATIC TRIANGLE
+      if (tags[evtri[0]]) {
+        meshDS->SetNodeOnFace(nodes[evtri[0]], TopoDS::Face(fmap(tag)));
+        tags[evtri[0]] = false;
       }
-      else {
-        tri = helper.AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
+      if (tags[evtri[1]]) {
+        meshDS->SetNodeOnFace(nodes[evtri[1]], TopoDS::Face(fmap(tag)));
+        tags[evtri[1]] = false;
       }
-      meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
+      if (tags[evtri[2]]) {
+        meshDS->SetNodeOnFace(nodes[evtri[2]], TopoDS::Face(fmap(tag)));
+        tags[evtri[2]] = false;
+      }
+      tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]],
+                            nodes[evtri[0]], nodes[evtri[1]], nodes[evtri[2]]);
+    }
+    else {
+      tri = helper.AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
     }
+    meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
+  }
 
-    /* enumerate quadrangles */
-    for(int it=1;it<=nq;it++) {
-      SMDS_MeshFace* quad;
-      mesh_get_quadrangle_vertices(msh, it, vtx);
-      mesh_get_quadrangle_extra_vertices(msh, it, &type, evquad);
-      mesh_get_quadrangle_tag(msh, it, &tag);
-      if (tags[vtx[0]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
-        tags[vtx[0]] = false;
-      };
-      if (tags[vtx[1]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
-        tags[vtx[1]] = false;
-      };
-      if (tags[vtx[2]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
-        tags[vtx[2]] = false;
-      };
-      if (tags[vtx[3]]) {
-        meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
-        tags[vtx[3]] = false;
-      };
-      if (type == MESHGEMS_MESH_ELEMENT_TYPE_QUAD9) {
-        // QUADRATIC QUADRANGLE
-        std::cout << "This is a quadratic quadrangle" << std::endl;
-        if (tags[evquad[0]]) {
-          meshDS->SetNodeOnFace(nodes[evquad[0]], TopoDS::Face(fmap(tag)));
-          tags[evquad[0]] = false;
-        }
-        if (tags[evquad[1]]) {
-          meshDS->SetNodeOnFace(nodes[evquad[1]], TopoDS::Face(fmap(tag)));
-          tags[evquad[1]] = false;
-        }
-        if (tags[evquad[2]]) {
-          meshDS->SetNodeOnFace(nodes[evquad[2]], TopoDS::Face(fmap(tag)));
-          tags[evquad[2]] = false;
-        }
-        if (tags[evquad[3]]) {
-          meshDS->SetNodeOnFace(nodes[evquad[3]], TopoDS::Face(fmap(tag)));
-          tags[evquad[3]] = false;
-        }
-        if (tags[evquad[4]]) {
-          meshDS->SetNodeOnFace(nodes[evquad[4]], TopoDS::Face(fmap(tag)));
-          tags[evquad[4]] = false;
-        }
-        quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]],
-                              nodes[evquad[0]], nodes[evquad[1]], nodes[evquad[2]], nodes[evquad[3]],
-                              nodes[evquad[4]]);
+  /* enumerate quadrangles */
+  for(int it=1;it<=nq;it++) {
+    SMDS_MeshFace* quad;
+    mesh_get_quadrangle_vertices(msh, it, vtx);
+    mesh_get_quadrangle_extra_vertices(msh, it, &type, evquad);
+    mesh_get_quadrangle_tag(msh, it, &tag);
+    if (tags[vtx[0]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
+      tags[vtx[0]] = false;
+    };
+    if (tags[vtx[1]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
+      tags[vtx[1]] = false;
+    };
+    if (tags[vtx[2]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
+      tags[vtx[2]] = false;
+    };
+    if (tags[vtx[3]]) {
+      meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
+      tags[vtx[3]] = false;
+    };
+    if (type == MESHGEMS_MESH_ELEMENT_TYPE_QUAD9) {
+      // QUADRATIC QUADRANGLE
+      std::cout << "This is a quadratic quadrangle" << std::endl;
+      if (tags[evquad[0]]) {
+        meshDS->SetNodeOnFace(nodes[evquad[0]], TopoDS::Face(fmap(tag)));
+        tags[evquad[0]] = false;
+      }
+      if (tags[evquad[1]]) {
+        meshDS->SetNodeOnFace(nodes[evquad[1]], TopoDS::Face(fmap(tag)));
+        tags[evquad[1]] = false;
       }
-      else {
-        quad = helper.AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
+      if (tags[evquad[2]]) {
+        meshDS->SetNodeOnFace(nodes[evquad[2]], TopoDS::Face(fmap(tag)));
+        tags[evquad[2]] = false;
       }
-      meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
+      if (tags[evquad[3]]) {
+        meshDS->SetNodeOnFace(nodes[evquad[3]], TopoDS::Face(fmap(tag)));
+        tags[evquad[3]] = false;
+      }
+      if (tags[evquad[4]]) {
+        meshDS->SetNodeOnFace(nodes[evquad[4]], TopoDS::Face(fmap(tag)));
+        tags[evquad[4]] = false;
+      }
+      quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]],
+                             nodes[evquad[0]], nodes[evquad[1]], nodes[evquad[2]], nodes[evquad[3]],
+                             nodes[evquad[4]]);
+    }
+    else {
+      quad = helper.AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
     }
+    meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
+  }
 
-    /* release the mesh object, the rest is released by cleaner */
-    cadsurf_data_regain_mesh(css, msh);
+  /* release the mesh object, the rest is released by cleaner */
+  cadsurf_data_regain_mesh(css, msh);
 
-    delete [] nodes;
-    delete [] tags;
+  delete [] nodes;
+  delete [] tags;
 
-    if ( needMerge ) // sew mesh computed by BLSURF with pre-existing mesh
+  if ( needMerge ) // sew mesh computed by BLSURF with pre-existing mesh
+  {
+    SMESH_MeshEditor editor( &aMesh );
+    SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
+    TIDSortedElemSet segementsOnEdge;
+    TIDSortedNodeSet nodesOnEdge;
+    TSubMeshSet::iterator smIt;
+    SMESHDS_SubMesh* smDS;
+    typedef SMDS_StdIterator< const SMDS_MeshNode*, SMDS_NodeIteratorPtr > TNodeIterator;
+    double tol;
+
+    // merge nodes on EDGE's with ones computed by BLSURF
+    for ( smIt = mergeSubmeshes.begin(); smIt != mergeSubmeshes.end(); ++smIt )
     {
-      SMESH_MeshEditor editor( &aMesh );
-      SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
-      TIDSortedElemSet segementsOnEdge;
-      TIDSortedNodeSet nodesOnEdge;
-      set< SMESHDS_SubMesh* >::iterator smIt;
-      SMESHDS_SubMesh* smDS;
-      typedef SMDS_StdIterator< const SMDS_MeshNode*, SMDS_NodeIteratorPtr > TNodeIterator;
-      double tol;
-
-      if ( !proxySubmeshes.empty() )
-      {
-        // merge proxy nodes with ones computed by BLSURF
-        for ( smIt = proxySubmeshes.begin(); smIt != proxySubmeshes.end(); ++smIt )
-        {
-          if (! (smDS = *smIt) ) continue;
-          nodesOnEdge.clear();
-          nodesOnEdge.insert( TNodeIterator( smDS->GetNodes() ), TNodeIterator() );
-          if (( smDS = meshDS->MeshElements( smDS->GetID() )))
-          {
-            nodesOnEdge.insert( TNodeIterator( smDS->GetNodes() ), TNodeIterator() );
-            if ( nodesOnEdge.size() > 1 )
-            {
-              tol = SMESH_TNodeXYZ( *nodesOnEdge.begin() ).Distance( *nodesOnEdge.rbegin() );
-              tol /= 50. * nodesOnEdge.size();
-              editor.FindCoincidentNodes( nodesOnEdge, tol, nodeGroupsToMerge );
-              editor.MergeNodes( nodeGroupsToMerge );
-            }
-          }
-        }
-      }
-      // merge nodes on EDGE's with ones computed by BLSURF
-      for ( smIt = mergeSubmeshes.begin(); smIt != mergeSubmeshes.end(); ++smIt )
-      {
-        if (! (smDS = *smIt) ) continue;
-        getNodeGroupsToMerge( smDS, meshDS->IndexToShape((*smIt)->GetID()), nodeGroupsToMerge );
+      if (! (smDS = *smIt) ) continue;
+      getNodeGroupsToMerge( smDS, meshDS->IndexToShape((*smIt)->GetID()), nodeGroupsToMerge );
 
-        SMDS_ElemIteratorPtr segIt = smDS->GetElements();
-        while ( segIt->more() )
-          segementsOnEdge.insert( segIt->next() );
-      }
-      // merge nodes
-      editor.MergeNodes( nodeGroupsToMerge );
+      SMDS_ElemIteratorPtr segIt = smDS->GetElements();
+      while ( segIt->more() )
+        segementsOnEdge.insert( segIt->next() );
+    }
+    // merge nodes
+    editor.MergeNodes( nodeGroupsToMerge );
 
-      // merge segments
-      SMESH_MeshEditor::TListOfListOfElementsID equalSegments;
-      editor.FindEqualElements( segementsOnEdge, equalSegments );
-      editor.MergeElements( equalSegments );
+    // merge segments
+    SMESH_MeshEditor::TListOfListOfElementsID equalSegments;
+    editor.FindEqualElements( segementsOnEdge, equalSegments );
+    editor.MergeElements( equalSegments );
 
-      // remove excess segments created on the boundary of viscous layers
-      const SMDS_TypeOfPosition onFace = SMDS_TOP_FACE;
-      for ( int i = 1; i <= emap.Extent(); ++i )
+    // remove excess segments created on the boundary of viscous layers
+    const SMDS_TypeOfPosition onFace = SMDS_TOP_FACE;
+    for ( int i = 1; i <= emap.Extent(); ++i )
+    {
+      if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( emap( i )))
       {
-        if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( emap( i )))
+        SMDS_ElemIteratorPtr segIt = smDS->GetElements();
+        while ( segIt->more() )
         {
-          SMDS_ElemIteratorPtr segIt = smDS->GetElements();
-          while ( segIt->more() )
-          {
-            const SMDS_MeshElement* seg = segIt->next();
-            if ( seg->GetNode(0)->GetPosition()->GetTypeOfPosition() == onFace ||
-                 seg->GetNode(1)->GetPosition()->GetTypeOfPosition() == onFace )
-              meshDS->RemoveFreeElement( seg, smDS );
-          }
+          const SMDS_MeshElement* seg = segIt->next();
+          if ( seg->GetNode(0)->GetPosition()->GetTypeOfPosition() == onFace ||
+               seg->GetNode(1)->GetPosition()->GetTypeOfPosition() == onFace )
+            meshDS->RemoveFreeElement( seg, smDS );
         }
       }
     }
-
-    if ( iLoop == 0 && !_haveViscousLayers)
-      break;
-
-  } // two loops
+  }
 
   // SetIsAlwaysComputed( true ) to sub-meshes of EDGEs w/o mesh
   TopLoc_Location loc; double f,l;
@@ -2403,6 +2461,15 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
       if ( !sm->IsMeshComputed() )
         sm->SetIsAlwaysComputed( true );
 
+  // Set error to FACE's w/o elements
+  for ( int i = 1; i <= fmap.Extent(); ++i )
+  {
+    SMESH_subMesh* sm = aMesh.GetSubMesh( fmap(i) );
+    if ( !sm->GetSubMeshDS() || sm->GetSubMeshDS()->NbElements() == 0 )
+      sm->GetComputeError().reset
+        ( new SMESH_ComputeError( COMPERR_ALGO_FAILED, _comment, this ));
+  }
+
   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 #ifndef WNT
   if ( oldFEFlags > 0 )
index 8ac0a45b859c3b827d7aaffbf54dcfe1aec23326..cbee8aeb2daf32b00810418d9a75ba6135faa667 100644 (file)
 #endif
 
 #include <Python.h>
-#include "SMESH_Algo.hxx"
-#include "SMESH_Mesh.hxx"
+#include <SMESH_Algo.hxx>
+#include <SMESH_Mesh.hxx>
 #include <SMESHDS_Mesh.hxx>
 #include <SMDS_MeshElement.hxx>
 #include <SMDS_MeshNode.hxx>
 #include <SMESH_Gen_i.hxx>
+#include <StdMeshers_ViscousLayers2D.hxx>
 #include <SALOMEconfig.h>
 #include CORBA_CLIENT_HEADER(SALOMEDS)
 #include CORBA_CLIENT_HEADER(GEOM_Gen)
@@ -89,7 +90,7 @@ class BLSURFPlugin_BLSURF: public SMESH_2D_Algo {
 
 #ifdef WITH_SMESH_CANCEL_COMPUTE
     virtual void CancelCompute();
-    bool computeCanceled() { return _compute_canceled;};
+    bool computeCanceled() { return _compute_canceled; }
 #endif
 
     virtual bool Evaluate(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape,
@@ -100,6 +101,9 @@ class BLSURFPlugin_BLSURF: public SMESH_2D_Algo {
     bool                           _haveViscousLayers;
 
   private:
+    bool compute(SMESH_Mesh&          aMesh,
+                 const TopoDS_Shape&  aShape);
+
     TopoDS_Shape entryToShape(std::string entry);
     void createEnforcedVertexOnFace(TopoDS_Shape FaceShape, BLSURFPlugin_Hypothesis::TEnfVertexList enfVertexList);
     void Set_NodeOnEdge(SMESHDS_Mesh* meshDS, SMDS_MeshNode* node, const TopoDS_Shape& ed);
index 1f3e4702ca4ab7701e647703d9b2ef2d22e30226..a8693fe8dffa12211568f6fb3e8a1f3a7c251c51 100644 (file)
 #ifndef BLSURFPLUGINGUI_H
 #define BLSURFPLUGINGUI_H
 
-#ifdef WIN32
-  #if defined BLSURFPLUGIN_GUI_DLG_EXPORTS || defined BLSURFPluginGUI_Dlg_EXPORTS
-    #define BLSURFPLUGIN_GUI_DLG_EXPORT __declspec( dllexport )
-  #else
-    #define BLSURFPLUGIN_GUI_DLG_EXPORT __declspec( dllimport )
-  #endif
-#else
-  #define BLSURFPLUGIN_GUI_DLG_EXPORT
-#endif
-
 enum PhysicalMesh
   {
     DefaultSize = 0,
@@ -67,8 +57,9 @@ enum {
 //////////////////////////////////////////
 
 #include "ui_BLSURFPluginGUI_StdWidget_QTD.h"
+#include "BLSURFPluginGUI_HypothesisCreator.h"
 
-class BLSURFPLUGIN_GUI_DLG_EXPORT BLSURFPluginGUI_StdWidget : public QWidget, 
+class BLSURFPLUGIN_GUI_EXPORT BLSURFPluginGUI_StdWidget : public QWidget, 
                                             public Ui::BLSURFPluginGUI_StdWidget_QTD
 {
   Q_OBJECT
@@ -90,7 +81,7 @@ public slots:
 
 #include "ui_BLSURFPluginGUI_AdvWidget_QTD.h"
 
-class BLSURFPLUGIN_GUI_DLG_EXPORT BLSURFPluginGUI_AdvWidget : public QWidget, 
+class BLSURFPLUGIN_GUI_EXPORT BLSURFPluginGUI_AdvWidget : public QWidget, 
                                             public Ui::BLSURFPluginGUI_AdvWidget_QTD
 {
   Q_OBJECT
index bb88c8969ea84e1baa0f8aba9c63685729c55132..e77ad36eff7523f4fdf09236f8c118a7810ee851 100644 (file)
@@ -73,9 +73,10 @@ libBLSURFPluginGUI_la_CPPFLAGS =     \
 libBLSURFPluginGUI_la_LDFLAGS =                        \
        $(QT_LIBS) $(QT_MT_LIBS)\
        ../BLSURFPlugin/libBLSURFEngine.la      \
+       ../../idl/libSalomeIDLBLSURFPLUGIN.la \
        $(GUI_LDFLAGS) -lqtx -lSalomeApp -lsuit -lSalomeObject -lLightApp       \
        $(GEOM_LDFLAGS) -lGEOM \
-       ${SMESH_LDFLAGS} -lSMESH -lGeomSelectionTools -lStdMeshersGUI -lSMESHFiltersSelection \
+       ${SMESH_LDFLAGS} -lSMESH  -lSMESHEngine -lGeomSelectionTools -lStdMeshersGUI -lSMESHFiltersSelection \
        $(CAS_KERNEL) $(MESHGEMS_CADSURF_LIBS)
 
 # resources files