Salome HOME
bos #16292 [CEA 656] MGCADSurf: option SetEnforced1D mesh
[plugins/blsurfplugin.git] / src / BLSURFPlugin / BLSURFPlugin_BLSURF.cxx
index 8f496d1a9018746319a8599228c600a680bd5b4c..e47c490c954f1095aed332128c14e4ed8f62c1b2 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2020  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2021  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // ---
 
 #include "BLSURFPlugin_BLSURF.hxx"
-#include "BLSURFPlugin_Hypothesis.hxx"
+
 #include "BLSURFPlugin_Attractor.hxx"
+#include "BLSURFPlugin_EnforcedMesh1D.hxx"
+#include "BLSURFPlugin_Hypothesis.hxx"
 
 extern "C"{
 #include <meshgems/meshgems.h>
 #include <meshgems/cadsurf.h>
 }
 
-#include <structmember.h>
-
-
-#include <Basics_Utils.hxx>
-
 #include <SMDS_EdgePosition.hxx>
 #include <SMESHDS_Group.hxx>
+#include <SMESH_File.hxx>
 #include <SMESH_Gen.hxx>
 #include <SMESH_Group.hxx>
 #include <SMESH_Mesh.hxx>
@@ -47,17 +45,20 @@ extern "C"{
 #include <SMESH_MesherHelper.hxx>
 #include <StdMeshers_FaceSide.hxx>
 #include <StdMeshers_ViscousLayers2D.hxx>
-#include <SMESH_File.hxx>
 
+#include <Basics_Utils.hxx>
 #include <utilities.h>
 
+#include <cstdlib>
 #include <limits>
 #include <list>
-#include <vector>
 #include <set>
-#include <cstdlib>
+#include <vector>
+
+#include <structmember.h> // python
 
 // OPENCASCADE includes
+#include <BRepBndLib.hxx>
 #include <BRepBuilderAPI_MakeFace.hxx>
 #include <BRepBuilderAPI_MakePolygon.hxx>
 #include <BRepBuilderAPI_MakeWire.hxx>
@@ -133,13 +134,13 @@ namespace
   static PyMethodDef PyStdOut_methods[] = {
     {"write",  (PyCFunction)PyStdOut_write,  METH_VARARGS,
       PyDoc_STR("write(string) -> None")},
-    {NULL,    NULL}   /* sentinel */
+    {0, 0, 0, 0}   /* sentinel */
   };
 
   static PyMemberDef PyStdOut_memberlist[] = {
     {(char*)"softspace", T_INT,  offsetof(PyStdOut, softspace), 0,
      (char*)"flag indicating that a space needs to be printed; used by print"},
-    {NULL} /* Sentinel */
+    {0, 0, 0, 0, 0} /* Sentinel */
   };
 
   static PyTypeObject PyStdOut_Type = {
@@ -187,6 +188,14 @@ namespace
     0,                            /*tp_new*/
     0,                            /*tp_free*/
     0,                            /*tp_is_gc*/
+    0,                            /*tp_bases*/
+    0,                            /*tp_mro*/
+    0,                            /*tp_cache*/
+    0,                            /*tp_subclasses*/
+    0,                            /*tp_weaklist*/
+    0,                            /*tp_del*/
+    0,                            /*tp_version_tag*/
+    0,                            /*tp_finalize*/
   };
 
   PyObject * newPyStdOut( std::string& out )
@@ -396,14 +405,12 @@ status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data);
 status_t size_on_vertex(integer vertex_id, real *size, void *user_data);
 
-typedef struct {
-        gp_XY uv;
-        gp_XYZ xyz;
-} projectionPoint;
-
 /////////////////////////////////////////////////////////
 
-projectionPoint getProjectionPoint(TopoDS_Face& theFace, const gp_Pnt& thePoint)
+BLSURFPlugin_BLSURF::projectionPoint
+BLSURFPlugin_BLSURF::getProjectionPoint(TopoDS_Face&  theFace,
+                                        const gp_Pnt& thePoint,
+                                        const bool    theAllowStateON)
 {
   projectionPoint myPoint;
 
@@ -431,16 +438,18 @@ projectionPoint getProjectionPoint(TopoDS_Face& theFace, const gp_Pnt& thePoint)
       {
         // check location on the face
         BRepClass_FaceClassifier FC( face, uv, BRep_Tool::Tolerance( face ));
-        if ( FC.State() == TopAbs_IN )
+        if (( FC.State() == TopAbs_IN ) ||
+            ( FC.State() == TopAbs_ON && theAllowStateON ))
         {
-          if ( !foundFace.IsNull() )
+          if ( !foundFace.IsNull() && !theAllowStateON )
             return myPoint; // thePoint seems to be TopAbs_ON
-          foundFace   = face;
-          myPoint.uv  = uv.XY();
-          myPoint.xyz = surface->Value( uv ).XYZ();
+          foundFace     = face;
+          myPoint.uv    = uv.XY();
+          myPoint.xyz   = surface->Value( uv ).XYZ();
+          myPoint.state = FC.State();
           // break;
         }
-        if ( FC.State() == TopAbs_ON )
+        if ( FC.State() == TopAbs_ON && !theAllowStateON )
           return myPoint;
       }
     }
@@ -453,12 +462,15 @@ projectionPoint getProjectionPoint(TopoDS_Face& theFace, const gp_Pnt& thePoint)
         const TopoDS_Face& face = d2f->second.first;
         const gp_Pnt2d &     uv = d2f->second.second;
         BRepClass_FaceClassifier FC( face, uv, Precision::Confusion());
-        if ( FC.State() == TopAbs_IN )
+        if (( FC.State() == TopAbs_IN ) ||
+            ( FC.State() == TopAbs_ON && theAllowStateON ))
         {
-          foundFace   = face;
-          myPoint.uv  = uv.XY();
-          myPoint.xyz = theHelper->GetSurface( face )->Value( uv ).XYZ();
-          break;
+          foundFace     = face;
+          myPoint.uv    = uv.XY();
+          myPoint.xyz   = theHelper->GetSurface( face )->Value( uv ).XYZ();
+          myPoint.state = FC.State();
+          if ( FC.State() == TopAbs_IN )
+            break;
         }
       }
     }
@@ -468,7 +480,7 @@ projectionPoint getProjectionPoint(TopoDS_Face& theFace, const gp_Pnt& thePoint)
     //                            "getProjectionPoint: can't find a face by a vertex");
     theFace = TopoDS::Face( foundFace );
   }
-  else
+  else // !theFace.IsNull()
   {
     Handle(Geom_Surface) surface = BRep_Tool::Surface( theFace );
     GeomAPI_ProjectPointOnSurf projector( thePoint, surface );
@@ -478,12 +490,14 @@ projectionPoint getProjectionPoint(TopoDS_Face& theFace, const gp_Pnt& thePoint)
 
     Standard_Real u,v;
     projector.LowerDistanceParameters(u,v);
-    myPoint.uv = gp_XY(u,v);
-    gp_Pnt aPnt = projector.NearestPoint();
-    myPoint.xyz = gp_XYZ(aPnt.X(),aPnt.Y(),aPnt.Z());
+    myPoint.uv  = gp_XY(u,v);
+    myPoint.xyz = projector.NearestPoint().XYZ();
 
     BRepClass_FaceClassifier FC( theFace, myPoint.uv, Precision::Confusion());
-    if ( FC.State() != TopAbs_IN )
+    myPoint.state = FC.State();
+
+    if (( FC.State() != TopAbs_IN ) &&
+        ( FC.State() != TopAbs_ON || !theAllowStateON ))
       theFace.Nullify();
   }
 
@@ -493,25 +507,20 @@ projectionPoint getProjectionPoint(TopoDS_Face& theFace, const gp_Pnt& thePoint)
 /////////////////////////////////////////////////////////
 TopoDS_Shape BLSURFPlugin_BLSURF::entryToShape(std::string entry)
 {
-  GEOM::GEOM_Object_var aGeomObj;
-  TopoDS_Shape S = TopoDS_Shape();
-  SALOMEDS::SObject_var aSObj = SMESH_Gen_i::getStudyServant()->FindObjectID( entry.c_str() );
-  if (!aSObj->_is_nil()) {
-    CORBA::Object_var obj = aSObj->GetObject();
-    aGeomObj = GEOM::GEOM_Object::_narrow(obj);
-    aSObj->UnRegister();
-  }
-  if ( !aGeomObj->_is_nil() )
-    S = SMESH_Gen_i::GetSMESHGen()->GeomObjectToShape( aGeomObj.in() );
+  GEOM::GEOM_Object_var aGeomObj = SMESH_Gen_i::GetSMESHGen()->GetGeomObjectByEntry( entry );
+  TopoDS_Shape                 S = SMESH_Gen_i::GetSMESHGen()->GeomObjectToShape( aGeomObj );
   return S;
 }
 
-void _createEnforcedVertexOnFace(TopoDS_Face faceShape, gp_Pnt aPnt, BLSURFPlugin_Hypothesis::TEnfVertex *enfVertex)
+void _createEnforcedVertexOnFace(TopoDS_Face                          faceShape,
+                                 const gp_Pnt&                        aPnt,
+                                 BLSURFPlugin_Hypothesis::TEnfVertex *enfVertex)
 {
   BLSURFPlugin_Hypothesis::TEnfVertexCoords enf_coords, coords, s_coords;
 
   // Find the face and get the (u,v) values of the enforced vertex on the face
-  projectionPoint myPoint = getProjectionPoint(faceShape,aPnt);
+  BLSURFPlugin_BLSURF::projectionPoint projPnt =
+    BLSURFPlugin_BLSURF::getProjectionPoint( faceShape, aPnt, /*allowStateON=*/true );
   if ( faceShape.IsNull() )
     return;
 
@@ -519,15 +528,16 @@ void _createEnforcedVertexOnFace(TopoDS_Face faceShape, gp_Pnt aPnt, BLSURFPlugi
   enf_coords.push_back(aPnt.Y());
   enf_coords.push_back(aPnt.Z());
 
-  coords.push_back(myPoint.uv.X());
-  coords.push_back(myPoint.uv.Y());
-  coords.push_back(myPoint.xyz.X());
-  coords.push_back(myPoint.xyz.Y());
-  coords.push_back(myPoint.xyz.Z());
+  coords.push_back(projPnt.uv.X());
+  coords.push_back(projPnt.uv.Y());
+  coords.push_back(projPnt.xyz.X());
+  coords.push_back(projPnt.xyz.Y());
+  coords.push_back(projPnt.xyz.Z());
+  coords.push_back(projPnt.state == TopAbs_ON);
 
-  s_coords.push_back(myPoint.xyz.X());
-  s_coords.push_back(myPoint.xyz.Y());
-  s_coords.push_back(myPoint.xyz.Z());
+  s_coords.push_back(projPnt.xyz.X());
+  s_coords.push_back(projPnt.xyz.Y());
+  s_coords.push_back(projPnt.xyz.Z());
 
   // Save pair projected vertex / enf vertex
   EnfVertexCoords2ProjVertex[s_coords] = enf_coords;
@@ -539,16 +549,10 @@ void _createEnforcedVertexOnFace(TopoDS_Face faceShape, gp_Pnt aPnt, BLSURFPlugi
     (*it)->grpName = enfVertex->grpName;
   }
 
-  int key = 0;
-  if (! FacesWithEnforcedVertices.Contains(faceShape)) {
-    key = FacesWithEnforcedVertices.Add(faceShape);
-  }
-  else {
-    key = FacesWithEnforcedVertices.FindIndex(faceShape);
-  }
+  int key = FacesWithEnforcedVertices.Add( faceShape );
 
   // If a node is already created by an attractor, do not create enforced vertex
-  int attractorKey = FacesWithSizeMap.FindIndex(faceShape);
+  int attractorKey = FacesWithSizeMap.FindIndex( faceShape );
   bool sameAttractor = false;
   if (attractorKey >= 0)
     if (FaceId2AttractorCoords.count(attractorKey) > 0)
@@ -557,7 +561,7 @@ void _createEnforcedVertexOnFace(TopoDS_Face faceShape, gp_Pnt aPnt, BLSURFPlugi
 
   if (FaceId2EnforcedVertexCoords.find(key) != FaceId2EnforcedVertexCoords.end()) {
     if (! sameAttractor)
-      FaceId2EnforcedVertexCoords[key].insert(coords); // there should be no redondant coords here (see std::set management)
+      FaceId2EnforcedVertexCoords[key].insert(coords); // there should be no redundant coords here (see std::set management)
   }
   else {
     if (! sameAttractor) {
@@ -567,17 +571,13 @@ void _createEnforcedVertexOnFace(TopoDS_Face faceShape, gp_Pnt aPnt, BLSURFPlugi
     }
   }
 }
-  
+
 /////////////////////////////////////////////////////////
 void BLSURFPlugin_BLSURF::createEnforcedVertexOnFace(TopoDS_Shape faceShape, BLSURFPlugin_Hypothesis::TEnfVertexList enfVertexList)
 {
-  BLSURFPlugin_Hypothesis::TEnfVertex* enfVertex;
   gp_Pnt aPnt;
-
-  BLSURFPlugin_Hypothesis::TEnfVertexList::const_iterator enfVertexListIt = enfVertexList.begin();
-
-  for( ; enfVertexListIt != enfVertexList.end() ; ++enfVertexListIt ) {
-    enfVertex = *enfVertexListIt;
+  for( BLSURFPlugin_Hypothesis::TEnfVertex* enfVertex : enfVertexList )
+  {
     // Case of manual coords
     if (enfVertex->coords.size() != 0) {
       aPnt.SetCoord(enfVertex->coords[0],enfVertex->coords[1],enfVertex->coords[2]);
@@ -586,7 +586,7 @@ void BLSURFPlugin_BLSURF::createEnforcedVertexOnFace(TopoDS_Shape faceShape, BLS
 
     // Case of geom vertex coords
     if (enfVertex->geomEntry != "") {
-      TopoDS_Shape GeomShape = entryToShape(enfVertex->geomEntry);
+      TopoDS_Shape     GeomShape = entryToShape(enfVertex->geomEntry);
       TopAbs_ShapeEnum GeomType  = GeomShape.ShapeType();
        if (GeomType == TopAbs_VERTEX)
        {
@@ -612,7 +612,7 @@ void BLSURFPlugin_BLSURF::createEnforcedVertexOnFace(TopoDS_Shape faceShape, BLS
 /////////////////////////////////////////////////////////
 void createAttractorOnFace(TopoDS_Shape GeomShape, std::string AttractorFunction, double defaultSize)
 {
-  double xa, ya, za; // Coordinates of attractor point
+  double xa=0., ya=0., za=0.; // Coordinates of attractor point
   double a, b;       // Attractor parameter
   double d = 0.;
   bool createNode=false; // To create a node on attractor projection
@@ -630,30 +630,30 @@ void createAttractorOnFace(TopoDS_Shape GeomShape, std::string AttractorFunction
   // xa
   pos1 = AttractorFunction.find(sep);
   if (pos1!=string::npos)
-  xa = atof(AttractorFunction.substr(10, pos1-10).c_str());
+    xa = atof(AttractorFunction.substr(10, pos1-10).c_str());
   // ya
   pos2 = AttractorFunction.find(sep, pos1+1);
   if (pos2!=string::npos) {
-  ya = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
-  pos1 = pos2;
+    ya = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
+    pos1 = pos2;
   }
   // za
   pos2 = AttractorFunction.find(sep, pos1+1);
   if (pos2!=string::npos) {
-  za = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
-  pos1 = pos2;
+    za = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
+    pos1 = pos2;
   }
   // a
   pos2 = AttractorFunction.find(sep, pos1+1);
   if (pos2!=string::npos) {
-  a = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
-  pos1 = pos2;
+    a = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
+    pos1 = pos2;
   }
   // b
   pos2 = AttractorFunction.find(sep, pos1+1);
   if (pos2!=string::npos) {
-  b = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
-  pos1 = pos2;
+    b = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
+    pos1 = pos2;
   }
   // createNode
   pos2 = AttractorFunction.find(sep, pos1+1);
@@ -665,13 +665,14 @@ void createAttractorOnFace(TopoDS_Shape GeomShape, std::string AttractorFunction
   // d
   pos2 = AttractorFunction.find(")");
   if (pos2!=string::npos) {
-  d = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
+    d = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
   }
 
   // Get the (u,v) values of the attractor on the face
-  projectionPoint myPoint = getProjectionPoint(TopoDS::Face(GeomShape),gp_Pnt(xa,ya,za));
-  gp_XY uvPoint = myPoint.uv;
-  gp_XYZ xyzPoint = myPoint.xyz;
+  BLSURFPlugin_BLSURF::projectionPoint projPnt =
+    BLSURFPlugin_BLSURF::getProjectionPoint(TopoDS::Face(GeomShape),gp_Pnt(xa,ya,za));
+  gp_XY    uvPoint = projPnt.uv;
+  gp_XYZ  xyzPoint = projPnt.xyz;
   Standard_Real u0 = uvPoint.X();
   Standard_Real v0 = uvPoint.Y();
   Standard_Real x0 = xyzPoint.X();
@@ -982,6 +983,14 @@ void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp,
       if ( !opIt->second.empty() ) {
         set_param(css, opIt->first.c_str(), opIt->second.c_str());
       }
+
+    if ( hyp->GetHyperPatches().size() < hyp->GetHyperPatchEntries().size() )
+    {
+      std::map< std::string, TopoDS_Shape > entryToShape;
+      FillEntryToShape( hyp, entryToShape );
+      const_cast<BLSURFPlugin_Hypothesis*>( hyp )->SetHyperPatchIDsByEntry( theGeomShape,
+                                                                            entryToShape );
+    }
   }
 
   if ( BLSURFPlugin_Hypothesis::HasPreCADOptions( hyp ))
@@ -994,7 +1003,7 @@ void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp,
   //set_param(css, "process_3d_topology",    _precadProcess3DTopology ? "1" : "0");
   //set_param(css, "discard_input_topology", _precadDiscardInput ? "1" : "0");
   //set_param(css, "max_number_of_points_per_patch", "1000000");
-  
+
    bool useGradation = false;
    switch (_physicalMesh)
    {
@@ -1468,7 +1477,7 @@ namespace
         //       distene_sizemap_delete(iso_sizemap_e);
         //     if(iso_sizemap_f)
         //       distene_sizemap_delete(iso_sizemap_f);
-        // 
+        //
         // //     if(clean_geo_sizemap_e)
         // //       distene_sizemap_delete(clean_geo_sizemap_e);
         // //     if(clean_geo_sizemap_f)
@@ -1548,7 +1557,7 @@ namespace
       //double tol = (( u2node.rbegin()->first - u2node.begin()->first ) / 20.) / u2node.size();
       Standard_Real f,l;
       BRep_Tool::Range( TopoDS::Edge( shape ), f,l );
-      double tol = (( l - f ) / 10.) / u2node.size(); // 10. - adjusted for #17262
+      double tol = (( l - f ) / 10.) / double( u2node.size() ); // 10. - adjusted for #17262
 
       std::multimap< double, const SMDS_MeshNode* >::iterator un2, un1;
       for ( un2 = u2node.begin(), un1 = un2++; un2 != u2node.end(); un1 = un2++ )
@@ -1602,7 +1611,7 @@ namespace
 
     TmpMesh()
     {
-      _myMeshDS = new SMESHDS_Mesh( _id, true );
+      _meshDS = new SMESHDS_Mesh( _id, true );
     }
     //--------------------------------------------------------------------------------
     /*!
@@ -1801,7 +1810,7 @@ namespace
           if ( !n2nIt->second ) {
             n->GetXYZ( xyz );
             gp_XY uv = tmpHelper.GetNodeUV( _proxyFace, n );
-            n2nIt->second = helper.AddNode( xyz[0], xyz[1], xyz[2], uv.X(), uv.Y() );
+            n2nIt->second = helper.AddNode( xyz[0], xyz[1], xyz[2], /*id=*/0, uv.X(), uv.Y() );
           }
           nodes[ nbN ] = n2nIt->second;
         }
@@ -1914,21 +1923,28 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
 
   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
   SMESH_MesherHelper helper( aMesh ), helperWithShape( aMesh );
-  myHelper = theHelper = & helperWithShape;
+  theHelper = & helperWithShape;
   // do not call helper.IsQuadraticSubMesh() because sub-meshes
   // may be cleaned and helper.myTLinkNodeMap gets invalid in such a case
   bool haveQuadraticSubMesh = helperWithShape.IsQuadraticSubMesh( aShape );
+  haveQuadraticSubMesh = haveQuadraticSubMesh || (_hypothesis && _hypothesis->GetQuadraticMesh());
+  helper.SetIsQuadratic( haveQuadraticSubMesh );
+  helperWithShape.SetIsQuadratic( haveQuadraticSubMesh );
+
   bool quadraticSubMeshAndViscousLayer = false;
   bool needMerge = false;
   typedef set< SMESHDS_SubMesh*, ShapeTypeCompare > TSubMeshSet;
   TSubMeshSet edgeSubmeshes;
   TSubMeshSet& mergeSubmeshes = edgeSubmeshes;
+  double existingPhySize = 0;
 
   TopTools_IndexedMapOfShape pmap, emap, fmap;
 
   TopTools_IndexedDataMapOfShapeListOfShape e2ffmap;
   TopExp::MapShapesAndAncestors( aShape, TopAbs_EDGE, TopAbs_FACE, e2ffmap );
 
+  BLSURFPlugin_EnforcedMesh1D enforcedMesh( helperWithShape, _hypothesis );
+
   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 #ifndef WIN32
   feclearexcept( FE_ALL_EXCEPT );
@@ -1970,9 +1986,6 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
 
   SetParameters(_hypothesis, css, aShape);
 
-  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 ) {
@@ -2012,8 +2025,8 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
   {
     TopoDS_Face f = TopoDS::Face(face_iter.Current());
 
-    SMESH_subMesh* fSM = aMesh.GetSubMesh( f );
-    if ( !fSM->IsEmpty() ) continue; // skip already meshed FACE with viscous layers
+    //SMESH_subMesh* fSM = aMesh.GetSubMesh( f );
+    //if ( !fSM->IsEmpty() ) continue; // skip already meshed FACE with viscous layers
 
     // make INTERNAL face oriented FORWARD (issue 0020993)
     if (f.Orientation() != TopAbs_FORWARD && f.Orientation() != TopAbs_REVERSED )
@@ -2070,35 +2083,29 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       }
 
       // 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)
+      std::map<int,std::vector<double> >::iterator attractor_iter =
+        FaceId2AttractorCoords.find( faceKey );
+      if (attractor_iter != FaceId2AttractorCoords.end() )
+      {
+        double * xyzCoords = & attractor_iter->second[2];
+        gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
+        BRepClass_FaceClassifier scl(f,P,1e-7);
+        TopAbs_State result = scl.State();
+        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 )
         {
-          double xyzCoords[3]  = {attractor_iter->second[2],
-                                  attractor_iter->second[3],
-                                  attractor_iter->second[4]};
-
-          gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
-          BRepClass_FaceClassifier scl(f,P,1e-7);
-          scl.Perform(f, P, 1e-7);
-          TopAbs_State result = scl.State();
-          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
-            double uvCoords[2] = {attractor_iter->second[0],attractor_iter->second[1]};
-            ienf++;
-            cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
-            cad_point_set_tag(point_p, ienf);
-          }
-          FaceId2AttractorCoords.erase(faceKey);
+          // Point is inside face and not on border
+          double * uvCoords = & attractor_iter->second[0];
+          ienf++;
+          cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
+          cad_point_set_tag(point_p, ienf);
         }
+        FaceId2AttractorCoords.erase(faceKey);
       }
 
       // -----------------
@@ -2128,14 +2135,19 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       BLSURFPlugin_Hypothesis::TEnfVertexCoordsList::const_iterator evlIt = evl.begin();
       for (; evlIt != evl.end(); ++evlIt)
       {
+        BLSURFPlugin_Hypothesis::TEnfVertexCoords xyzCoords = { evlIt->at(2),
+                                                                evlIt->at(3),
+                                                                evlIt->at(4)};
+        bool isOnEdge = evlIt->at(5);
+        if ( isOnEdge )
+        {
+          enforcedMesh.AddVertexOnEdge( xyzCoords.data() );
+          continue;
+        }
         double uvCoords[2] = { evlIt->at(0), evlIt->at(1) };
         ienf++;
         cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
         int tag = 0;
-        BLSURFPlugin_Hypothesis::TEnfVertexCoords xyzCoords;
-        xyzCoords.push_back(evlIt->at(2));
-        xyzCoords.push_back(evlIt->at(3));
-        xyzCoords.push_back(evlIt->at(4));
         std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(xyzCoords);
         if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end() &&
             !enfCoordsIt->second.empty() )
@@ -2164,19 +2176,24 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
                                            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())
+
+    std::vector< TopoDS_Edge > edges;
+    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 TopoDS_Edge& e = TopoDS::Edge( edge_iter.Current() );
+      if ( !enforcedMesh.GetSplitsOfEdge( e, edges, emap ))
+        edges.push_back( e );
+    }
+    for ( const TopoDS_Edge& e : edges )
+    {
+      int ic = emap.Add(e);
 
       double tmin,tmax;
       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
 
-      if (HasSizeMapOnEdge){
-        edgeKey = EdgesWithSizeMap.FindIndex(e);
+      if ( HasSizeMapOnEdge )
+      {
+        int edgeKey = EdgesWithSizeMap.FindIndex(e);
         if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end())
         {
           theSizeMapStr = EdgeId2SizeMap[edgeKey];
@@ -2221,6 +2238,8 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
           // tmin and tmax can change in case of viscous layer on an adjacent edge
           tmin = nodeDataVec.front().param;
           tmax = nodeDataVec.back().param;
+
+          existingPhySize += nodeData->Length() / double( nodeDataVec.size() - 1 );
         }
         else
         {
@@ -2236,6 +2255,7 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
          The following call sets it to the same value as the edge_id : */
       // IMP23368. Do not set tag to an EDGE shared by FACEs of a hyper-patch
       bool isInHyperPatch = false;
+      if ( e2ffmap.Contains( e )) // not there for a part of EDGE divided by enforced segments
       {
         std::set< int > faceTags, faceIDs;
         TopTools_ListIteratorOfListOfShape fIt( e2ffmap.FindFromKey( e ));
@@ -2255,7 +2275,7 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       if ( !isInHyperPatch )
         cad_edge_set_tag(edg, ic);
 
-      /* by default, an edge does not necessalry appear in the resulting mesh,
+      /* by default, an edge does not necessarily appear in the resulting mesh,
          unless the following property is set :
       */
       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
@@ -2268,7 +2288,7 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       if ( nodeData )
       {
         const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
-        const int                      nbNodes     = nodeDataVec.size();
+        const int                      nbNodes     = (int) nodeDataVec.size();
 
         dcad_edge_discretization_t *dedge;
         dcad_get_edge_discretization(dcad, edg, &dedge);
@@ -2279,9 +2299,9 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
         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 );
+          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
@@ -2304,38 +2324,23 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       *****************************************************************************************/
 
       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);
-          // SMESH_subMesh* sm = aMesh.GetSubMesh(v);
-          // if ( sm->IsMeshComputed() )
-          //   edgeSubmeshes.insert( sm->GetSubMeshDS() );
-        }
+      int   ip[2];
+      double d[2];
+      gp_Pnt2d uv0 = curves.back()->Value(tmin);
+      gp_Pnt    p0 = surfaces.back()->Value( uv0.X(), uv0.Y() );
 
-//        std::string aFileName = "fmap_vertex_";
-//        aFileName.append(val_to_string(*ip));
-//        aFileName.append(".brep");
-//        BRepTools::Write(v,aFileName.c_str());
+      for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next())
+      {
+        const TopoDS_Vertex& v = TopoDS::Vertex(ex_edge.Current());
+        ip[ npts ] = pmap.Add(v);
+        d [ npts ] = p0.SquareDistance(BRep_Tool::Pnt(v));
+        ++npts;
 
-        if (HasSizeMapOnVertex){
-          vertexKey = VerticesWithSizeMap.FindIndex(v);
-          if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
+        if (HasSizeMapOnVertex)
+        {
+          int vertexKey = VerticesWithSizeMap.FindIndex(v);
+          if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end())
+          {
             theSizeMapStr = VertexId2SizeMap[vertexKey];
             if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
               continue;
@@ -2351,62 +2356,66 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
             PyGILState_Release(gstate);
           }
         }
-      }
-      if (npts != 2) {
-        // should not happen
+      } // loop on vertices
+
+      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);
-        }
+        continue;
       }
+
+      if ( d[0] > d[1] )
+        std::swap( ip[0], ip[1] );
+
+      // This defines the curves extremity connectivity
+      cad_edge_set_extremities    (edg, ip[0], ip[1]);
+      /* set the tag (color) to the same value as the extremity id : */
+      cad_edge_set_extremities_tag(edg, ip[0], ip[1]);
+
     } // 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() )
+    // ==============================
+    // Add segments of enforced mesh
+    // ==============================
+
+    if ( enforcedMesh.HasSegmentsOnFace( f ))
     {
-      const SMDS_MeshNode* n = nIt->next();
-      if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
+      BLSURFPlugin_EnforcedMesh1D::Segmemnt seg;
+      while ( enforcedMesh.NextSegment( seg, pmap ))
       {
-        needMerge = true; // to correctly sew with viscous mesh
-        // 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()));
-        }
-        continue;
+        curves.push_back( seg._pcurve );
+
+        cad_edge_t *edg = cad_edge_new( fce, seg._tag, seg._u[0], seg._u[1],
+                                        curv_fun,seg._pcurve.get());
+        cad_edge_set_tag( edg, seg._tag );
+
+        cad_edge_set_property( edg, EDGE_PROPERTY_SOFT_REQUIRED );
+        cad_edge_set_property( edg, EDGE_PROPERTY_INTERNAL );
+
+        cad_edge_set_extremities    ( edg, seg._vTag[0], seg._vTag[1]);
+        cad_edge_set_extremities_tag( edg, seg._vTag[0], seg._vTag[1]);
+
+        dcad_edge_discretization_t *dedge;
+        dcad_get_edge_discretization( dcad, edg, &dedge );
+        dcad_edge_discretization_set_vertex_count( dedge, 2 );
+
+        dcad_edge_discretization_set_vertex_coordinates( dedge, 1,
+                                                         seg._u  [0],
+                                                         &seg._uv[0].ChangeCoord(1),
+                                                         seg._xyz[0].ChangeData() );
+        dcad_edge_discretization_set_vertex_coordinates( dedge, 2,
+                                                         seg._u  [1],
+                                                         &seg._uv[1].ChangeCoord(1),
+                                                         seg._xyz[1].ChangeData() );
+
+        dcad_edge_discretization_set_vertex_tag( dedge, 1, seg._vTag[0] );
+        dcad_edge_discretization_set_vertex_tag( dedge, 2, seg._vTag[1] );
+
+        dcad_edge_discretization_set_property(dedge, DISTENE_DCAD_PROPERTY_REQUIRED);
       }
     }
-    if ( allowSubMeshClearing )
-    {
-      SMDS_ElemIteratorPtr eIt = smDS->GetElements();
-      while ( eIt->more() ) meshDS->RemoveFreeElement( eIt->next(), 0 );
-      SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
-      while ( nIt->more() ) meshDS->RemoveFreeNode( nIt->next(), 0 );
-      smDS->Clear();
-    }
-    else
-    {
-      needMerge = true;
-    }
-  }
+
+
+  } //for face
 
   ///////////////////////
   // PERIODICITY       //
@@ -2433,8 +2442,9 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       {
         // If no source points, call periodicity without transformation function
         meshgems_cad_periodicity_transformation_t periodicity_transformation = NULL;
-        status = cad_add_face_multiple_periodicity_with_transformation_function(c, theFace1_ids_c, theFace1_ids.size(),
-                                                                                theFace2_ids_c, theFace2_ids.size(), periodicity_transformation, NULL);
+        status = cad_add_face_multiple_periodicity_with_transformation_function
+          (c, theFace1_ids_c, (meshgems_integer) theFace1_ids.size(), theFace2_ids_c,
+           (meshgems_integer) theFace2_ids.size(), periodicity_transformation, NULL);
         if(status != STATUS_OK)
           cout << "cad_add_face_multiple_periodicity_with_transformation_function failed with error code " << status << "\n";
       }
@@ -2443,11 +2453,12 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
         // get the transformation vertices
         double* theSourceVerticesCoords_c = &_preCadFacesIDsPeriodicityVector[i].theSourceVerticesCoords[0];
         double* theTargetVerticesCoords_c = &_preCadFacesIDsPeriodicityVector[i].theTargetVerticesCoords[0];
-        int nbSourceVertices = _preCadFacesIDsPeriodicityVector[i].theSourceVerticesCoords.size()/3;
-        int nbTargetVertices = _preCadFacesIDsPeriodicityVector[i].theTargetVerticesCoords.size()/3;
+        int nbSourceVertices = (int)_preCadFacesIDsPeriodicityVector[i].theSourceVerticesCoords.size()/3;
+        int nbTargetVertices = (int)_preCadFacesIDsPeriodicityVector[i].theTargetVerticesCoords.size()/3;
 
-        status = cad_add_face_multiple_periodicity_with_transformation_function_by_points(c, theFace1_ids_c, theFace1_ids.size(),
-                                                                                          theFace2_ids_c, theFace2_ids.size(), theSourceVerticesCoords_c, nbSourceVertices, theTargetVerticesCoords_c, nbTargetVertices);
+        status = cad_add_face_multiple_periodicity_with_transformation_function_by_points
+          (c, theFace1_ids_c, (meshgems_integer) theFace1_ids.size(), theFace2_ids_c,
+           (meshgems_integer) theFace2_ids.size(), theSourceVerticesCoords_c, nbSourceVertices, theTargetVerticesCoords_c, nbTargetVertices);
         if(status != STATUS_OK)
           cout << "cad_add_face_multiple_periodicity_with_transformation_function_by_points failed with error code " << status << "\n";
       }
@@ -2478,8 +2489,9 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
       {
         // If no source points, call periodicity without transformation function
         meshgems_cad_periodicity_transformation_t periodicity_transformation = NULL;
-        status = cad_add_edge_multiple_periodicity_with_transformation_function(c, theEdge1_ids_c, theEdge1_ids.size(),
-                                                                                theEdge2_ids_c, theEdge2_ids.size(), periodicity_transformation, NULL);
+        status = cad_add_edge_multiple_periodicity_with_transformation_function
+          (c, theEdge1_ids_c, (meshgems_integer) theEdge1_ids.size(), theEdge2_ids_c,
+           (meshgems_integer) theEdge2_ids.size(), periodicity_transformation, NULL);
         if(status != STATUS_OK)
           cout << "cad_add_edge_multiple_periodicity_with_transformation_function failed with error code " << status << "\n";
       }
@@ -2488,22 +2500,35 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
         // get the transformation vertices
         double* theSourceVerticesCoords_c = &_preCadEdgesIDsPeriodicityVector[i].theSourceVerticesCoords[0];
         double* theTargetVerticesCoords_c = &_preCadEdgesIDsPeriodicityVector[i].theTargetVerticesCoords[0];
-        int nbSourceVertices = _preCadEdgesIDsPeriodicityVector[i].theSourceVerticesCoords.size()/3;
-        int nbTargetVertices = _preCadEdgesIDsPeriodicityVector[i].theTargetVerticesCoords.size()/3;
+        int nbSourceVertices = (int) _preCadEdgesIDsPeriodicityVector[i].theSourceVerticesCoords.size()/3;
+        int nbTargetVertices = (int) _preCadEdgesIDsPeriodicityVector[i].theTargetVerticesCoords.size()/3;
 
-        status = cad_add_edge_multiple_periodicity_with_transformation_function_by_points(c, theEdge1_ids_c, theEdge1_ids.size(),
-                                                                                          theEdge2_ids_c, theEdge2_ids.size(), theSourceVerticesCoords_c, nbSourceVertices, theTargetVerticesCoords_c, nbTargetVertices);
+        status = cad_add_edge_multiple_periodicity_with_transformation_function_by_points
+          (c, theEdge1_ids_c, (meshgems_integer) theEdge1_ids.size(), theEdge2_ids_c,
+           (meshgems_integer) theEdge2_ids.size(), theSourceVerticesCoords_c, nbSourceVertices, theTargetVerticesCoords_c, nbTargetVertices);
         if(status != STATUS_OK)
           cout << "cad_add_edge_multiple_periodicity_with_transformation_function_by_points failed with error code " << status << "\n";
       }
     }
   }
 
-  
-  // TODO: be able to use a mesh in input.
-  // See imsh usage in Products/templates/mg-cadsurf_template_common.cpp
-  // => cadsurf_set_mesh
-    
+  if ( !_hypothesis && !edgeSubmeshes.empty() && existingPhySize != 0 )
+  {
+    // prevent failure due to the default PhySize incompatible with size of existing 1D mesh
+    // and with face size
+    // double minFaceSize = existingPhySize / edgeSubmeshes.size();
+    // for ( int iF = 1; iF <= fmap.Extent(); ++iF )
+    // {
+    //   Bnd_Box box;
+    //   BRepBndLib::Add( fmap( iF ), box );
+    //   gp_XYZ delta = box.CornerMax().XYZ() - box.CornerMin().XYZ();
+    //   std::sort( delta.ChangeData(), delta.ChangeData() + 3 );
+    //   minFaceSize = Min( minFaceSize, delta.Coord(2) );
+    // }
+    // set_param(css, "global_physical_size", val_to_string( minFaceSize * 0.5 ).c_str());
+    // set_param(css, "max_size",             val_to_string( minFaceSize * 5 ).c_str());
+  }
+
   // Use the original dcad
   cadsurf_set_dcad(css, dcad);
 
@@ -2541,16 +2566,56 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
 
   mesh_t *msh = NULL;
   cadsurf_get_mesh(css, &msh);
-  if(!msh){
+  if ( !msh || STATUS_IS_ERROR( status ))
+  {
     /* release the mesh object */
     cadsurf_regain_mesh(css, msh);
     return error(_comment);
   }
 
+  // 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; // to correctly sew with viscous mesh
+        // 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()));
+        }
+        continue;
+      }
+    }
+    if ( allowSubMeshClearing )
+    {
+      SMDS_ElemIteratorPtr eIt = smDS->GetElements();
+      while ( eIt->more() ) meshDS->RemoveFreeElement( eIt->next(), 0 );
+      SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
+      while ( nIt->more() ) meshDS->RemoveFreeNode( nIt->next(), 0 );
+      smDS->Clear();
+    }
+    else
+    {
+      needMerge = true;
+    }
+  }
+
   std::string GMFFileName = BLSURFPlugin_Hypothesis::GetDefaultGMFFile();
   if (_hypothesis)
     GMFFileName = _hypothesis->GetGMFFile();
-  if (GMFFileName != "") {
+  if (GMFFileName != "")
+  {
     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)
@@ -2563,9 +2628,9 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
   integer *evedg, *evtri, *evquad, *tags_buff, type;
   real xyz[3];
 
-  mesh_get_vertex_count(msh, &nv);
-  mesh_get_edge_count(msh, &ne);
-  mesh_get_triangle_count(msh, &nt);
+  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);
@@ -2576,82 +2641,80 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
   std::vector<const SMDS_MeshNode*> nodes(nv+1);
   std::vector<bool>                  tags(nv+1);
 
+  BLSURFPlugin_Hypothesis::TEnfVertexCoords projVertex;
+
   /* enumerated vertices */
-  for(int iv=1;iv<=nv;iv++) {
+  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
     nodes[iv] = NULL;
-    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]), 1e3*tol))
-        xyz[0] = p.X(), xyz[1] = p.Y(), xyz[2] = p.Z();
-      else
-        tag = 0; // enforced or attracted vertex
-      nodes[iv] = SMESH_Algo::VertexNode( v, meshDS );
+    bool isEnforcedNode = false;
+    if ( tag > 0 )
+    {
+      if ( 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]), 1e3*tol))
+          xyz[0] = p.X(), xyz[1] = p.Y(), xyz[2] = p.Z();
+        else
+          tag = 0; // enforced or attracted vertex
+        nodes[iv] = SMESH_Algo::VertexNode( v, meshDS );
+      }
+      if ( !nodes[iv] )
+        isEnforcedNode = ( nodes[iv] = enforcedMesh.GetNodeByTag( tag, pmap ));
     }
     if ( !nodes[iv] )
       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);
+    projVertex.assign( xyz, xyz + 3 );
+    auto enfCoordsIt = EnfVertexCoords2EnfVertexList.find( projVertex );
     if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end())
     {
-      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();
-          while (grIt->more()) {
+      for ( BLSURFPlugin_Hypothesis::TEnfVertex *currentEnfVertex : enfCoordsIt->second )
+        if (currentEnfVertex->grpName != "")
+        {
+          SMDS_MeshGroup* groupDS = nullptr;
+          for ( SMESH_Mesh::GroupIteratorPtr grIt = aMesh.GetGroups(); grIt->more() && !groupDS; )
+          {
             SMESH_Group * group = grIt->next();
-            if ( !group ) continue;
-            SMESHDS_GroupBase* groupDS = group->GetGroupDS();
-            if ( !groupDS ) continue;
-            if ( groupDS->GetType()==SMDSAbs_Node && currentEnfVertex->grpName.compare(group->GetName())==0) {
-              SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
-              aGroupDS->SMDSGroup().Add(nodes[iv]);
-              // How can I inform the hypothesis ?
-              //                 _hypothesis->AddEnfVertexNodeID(currentEnfVertex->grpName,nodes[iv]->GetID());
-              groupDone = true;
-              break;
-            }
+            SMESHDS_Group * grp = dynamic_cast< SMESHDS_Group* >( group->GetGroupDS() );
+            if ( grp &&
+                 grp->GetType() == SMDSAbs_Node &&
+                 currentEnfVertex->grpName == group->GetName() )
+              groupDS = &grp->SMDSGroup();
           }
-          if (!groupDone)
+          if ( !groupDS )
           {
-            SMESH_Group* aGroup = aMesh.AddGroup( SMDSAbs_Node, currentEnfVertex->grpName.c_str() );
-            aGroup->SetName( currentEnfVertex->grpName.c_str() );
-            SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
-            aGroupDS->SMDSGroup().Add(nodes[iv]);
-            groupDone = true;
+            SMESH_Group * group = aMesh.AddGroup( SMDSAbs_Node, currentEnfVertex->grpName.c_str() );
+            SMESHDS_Group * grp = static_cast< SMESHDS_Group* >( group->GetGroupDS() );
+            groupDS = &grp->SMDSGroup();
           }
-          if (!groupDone)
-            throw SALOME_Exception(LOCALIZED("An enforced vertex node was not added to a group"));
+          groupDS->Add( nodes[iv] );
         }
-        else
-          MESSAGE("Group name is empty: '"<<currentEnfVertex->grpName<<"' => group is not created");
-      }
     }
+    if ( isEnforcedNode )
+      continue;
 
     // internal points are tagged to zero
-    if(tag > 0 && tag <= pmap.Extent() ){
+    if ( tag > 0 && tag <= pmap.Extent() )
+    {
       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
       tags[iv] = false;
-    } else {
+    }
+    else
+    {
       tags[iv] = true;
     }
   }
 
   /* enumerate edges */
-  for(int it=1;it<=ne;it++) {
+  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);
@@ -2659,17 +2722,24 @@ bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&         aMesh,
 
     // If PreCAD performed some cleaning operations (remove tiny edges,
     // merge edges ...) an output tag can indeed represent several original tags.
-    // Get the initial tags corresponding to the output tag and redefine the tag as 
+    // Get the initial tags corresponding to the output tag and redefine the tag as
     // the last of the two initial tags (else the output tag is out of emap and hasn't any meaning)
     mesh_get_composite_tag_definition(msh, tag, &nb_tag, tags_buff);
-    if(nb_tag > 1)  
-      tag=tags_buff[nb_tag-1];
+    if ( nb_tag > 1 )
+      tag = tags_buff[ nb_tag-1 ];
     if ( tag < 1 || tag > emap.Extent() )
     {
-      std::cerr << "MG-CADSurf BUG:::: Edge tag " << tag
-                << " does not point to a CAD edge (nb edges " << emap.Extent() << ")" << std::endl;
+      if ( !enforcedMesh.IsSegmentTag( tag )) // it's a false INTERNAL EDGE of enforced mesh
+      {
+        std::cerr << "MG-CADSurf BUG:::: Edge tag " << tag
+                  << " does not point to a CAD edge (nb edges " << emap.Extent() << ")"
+                  << std::endl;
+      }
       continue;
     }
+    if ( meshDS->ShapeToIndex( emap( tag )) == 0 )
+      tag = enforcedMesh.GetTagOfSplitEdge( tag );
+
     if (tags[vtx[0]]) {
       Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
       tags[vtx[0]] = false;
@@ -2923,6 +2993,10 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelpe
   context_t *ctx = context_new();
   if (!ctx) return error("Pb in context_new()");
 
+  if ( aMesh.NbNodes() > std::numeric_limits< meshgems_integer >::max() ||
+       aMesh.NbFaces() > std::numeric_limits< meshgems_integer >::max() )
+    return error("Too large input mesh");
+
   BLSURF_Cleaner cleaner( ctx );
 
   message_cb_user_data mcud;
@@ -2941,7 +3015,7 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelpe
   // Fill an input mesh
 
   mesh_t * msh = meshgems_mesh_new_in_memory( ctx );
-  if ( !msh ) return error("Pb. in meshgems_mesh_new_in_memory()"); 
+  if ( !msh ) return error("Pb. in meshgems_mesh_new_in_memory()");
 
   // mark nodes used by 2D elements
   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
@@ -2951,7 +3025,7 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelpe
     const SMDS_MeshNode* n = nodeIt->next();
     n->setIsMarked( n->NbInverseElements( SMDSAbs_Face ));
   }
-  meshgems_mesh_set_vertex_count( msh, meshDS->NbNodes() );
+  meshgems_mesh_set_vertex_count( msh, (meshgems_integer) meshDS->NbNodes() );
 
   // set node coordinates
   if ( meshDS->NbNodes() != meshDS->MaxNodeID() )
@@ -2968,8 +3042,8 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelpe
   }
 
   // set nodes of faces
-  meshgems_mesh_set_triangle_count  ( msh, meshDS->GetMeshInfo().NbTriangles() );
-  meshgems_mesh_set_quadrangle_count( msh, meshDS->GetMeshInfo().NbQuadrangles() );
+  meshgems_mesh_set_triangle_count  ( msh, (meshgems_integer) meshDS->GetMeshInfo().NbTriangles() );
+  meshgems_mesh_set_quadrangle_count( msh, (meshgems_integer) meshDS->GetMeshInfo().NbQuadrangles() );
   meshgems_integer nodeIDs[4];
   meshgems_integer iT = 1, iQ = 1;
   SMDS_FaceIteratorPtr faceIt = meshDS->facesIterator();
@@ -2980,7 +3054,7 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelpe
     if ( nbNodes > 4 || face->IsPoly() ) continue;
 
     for ( i = 0; i < nbNodes; ++i )
-      nodeIDs[i] = face->GetNode( i )->GetID();
+      nodeIDs[i] = (meshgems_integer) face->GetNode( i )->GetID();
     if ( nbNodes == 3 )
       meshgems_mesh_set_triangle_vertices  ( msh, iT++, nodeIDs );
     else
@@ -3040,7 +3114,7 @@ bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelpe
   {
     meshgems_mesh_get_vertex_coordinates( omsh, i, xyz );
     SMDS_MeshNode* n = meshDS->AddNode( xyz[0], xyz[1], xyz[2] );
-    nodeID = n->GetID();
+    nodeID = (meshgems_integer) n->GetID();
     meshgems_mesh_set_vertex_tag( omsh, i, &nodeID ); // save mapping of IDs in MG and SALOME meshes
   }
 
@@ -3162,31 +3236,27 @@ void BLSURFPlugin_BLSURF::Set_NodeOnEdge(SMESHDS_Mesh*        meshDS,
  */
 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
 {
-  /* t is given. It contains the t (time) 1D parametric coordintaes
-     of the point PreCAD/MG-CADSurf is querying on the curve */
+  /* t is given. It contains the 1D parametric coordintaes
+     of the point MG-CADSurf is querying on the curve */
 
-  /* user_data identifies the edge PreCAD/MG-CADSurf is querying
-   * (see cad_edge_new later in this example) */
-  const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
+  /* user_data identifies the edge MG-CADSurf is querying */
+  const Geom2d_Curve* pargeo = (const Geom2d_Curve*) user_data;
 
-  if (uv){
-   /* MG-CADSurf is querying the function evaluation */
-    gp_Pnt2d P;
-    P=pargeo->Value(t);
+  if (uv) {
+    /* MG-CADSurf is querying the function evaluation */
+    gp_Pnt2d P = pargeo->Value(t);
     uv[0]=P.X(); uv[1]=P.Y();
   }
 
-  if(dt) {
-   /* query for the first order derivatives */
-    gp_Vec2d V1;
-    V1=pargeo->DN(t,1);
+  if (dt) {
+    /* query for the first order derivatives */
+    gp_Vec2d V1 = pargeo->DN(t,1);
     dt[0]=V1.X(); dt[1]=V1.Y();
   }
 
-  if(dtt){
+  if (dtt) {
     /* query for the second order derivatives */
-    gp_Vec2d V2;
-    V2=pargeo->DN(t,2);
+    gp_Vec2d V2 = pargeo->DN(t,2);
     dtt[0]=V2.X(); dtt[1]=V2.Y();
   }
 
@@ -3380,7 +3450,7 @@ status_t message_cb(message_t *msg, void *user_data)
        err.find("periodicity") != string::npos )
   {
     // remove ^A from the tail
-    int len = strlen( desc );
+    size_t len = strlen( desc );
     while (len > 0 && desc[len-1] != '\n')
       len--;
     mcud->_error->append( desc, len );
@@ -3486,8 +3556,8 @@ bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh&         aMesh,
       nb1d = (int)( fullAng/_angleMesh + 1 );
     }
     fullNbSeg += nb1d;
-    std::vector<int> aVec(SMDSEntity_Last);
-    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
+    std::vector<smIdType> aVec(SMDSEntity_Last);
+    for(smIdType i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
     if( IsQuadratic > 0 ) {
       aVec[SMDSEntity_Node] = 2*nb1d - 1;
       aVec[SMDSEntity_Quad_Edge] = nb1d;
@@ -3534,7 +3604,7 @@ bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh&         aMesh,
         nbTria = nbQuad = nbTria / 3 + 1;
       }
     }
-    std::vector<int> aVec(SMDSEntity_Last,0);
+    std::vector<smIdType> aVec(SMDSEntity_Last,0);
     if( IsQuadratic ) {
       int nb1d_in = (nbTria*3 - nb1d) / 2;
       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
@@ -3558,8 +3628,8 @@ bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh&         aMesh,
   double tetrVol = 0.1179*ELen*ELen*ELen;
   int nbVols  = int(aVolume/tetrVol);
   int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
-  std::vector<int> aVec(SMDSEntity_Last);
-  for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
+  std::vector<smIdType> aVec(SMDSEntity_Last);
+  for(smIdType i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
   if( IsQuadratic ) {
     aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
     aVec[SMDSEntity_Quad_Tetra] = nbVols;
@@ -3573,3 +3643,23 @@ bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh&         aMesh,
 
   return true;
 }
+
+//================================================================================
+/*!
+ * \brief Find TopoDS_Shape for each hyper-patch study entry in a hypothesis
+ */
+//================================================================================
+
+void BLSURFPlugin_BLSURF::FillEntryToShape( const BLSURFPlugin_Hypothesis*          hyp,
+                                            std::map< std::string, TopoDS_Shape > & entryToShape )
+{
+  SMESH_Gen_i* smeshGen = SMESH_Gen_i::GetSMESHGen();
+  for ( const ::BLSURFPlugin_Hypothesis::THyperPatchEntries& entries : hyp->GetHyperPatchEntries() )
+    for ( const std::string& entry : entries )
+    {
+      GEOM::GEOM_Object_var go = smeshGen->GetGeomObjectByEntry( entry );
+      TopoDS_Shape       shape = smeshGen->GeomObjectToShape( go );
+      if ( !shape.IsNull() )
+        entryToShape.insert({ entry, shape });
+    }
+}