src/BLSURFPlugin/BLSURFPlugin_BLSURF.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18
  19 // ---
  20 // File    : BLSURFPlugin_BLSURF.cxx
  21 // Authors : Francis KLOSS (OCC) & Patrick LAUG (INRIA) & Lioka RAZAFINDRAZAKA (CEA)
  22 //           & Aurelien ALLEAUME (DISTENE)
  23 //           Size maps developement: Nicolas GEIMER (OCC) & Gilles DAVID (EURIWARE)
  24 // ---
  25
  26 #include "BLSURFPlugin_BLSURF.hxx"
  27 #include "BLSURFPlugin_Hypothesis.hxx"
  28 #include "BLSURFPlugin_Attractor.hxx"
  29
  30 extern "C"{
  31 #include <meshgems/meshgems.h>
  32 #include <meshgems/cadsurf.h>
  33 #include <meshgems/precad.h>
  34 }
  35
  36 #include <structmember.h>
  37
  38
  39 #include <Basics_Utils.hxx>
  40 #include <Basics_OCCTVersion.hxx>
  41
  42 #include <SMDS_EdgePosition.hxx>
  43 #include <SMESHDS_Group.hxx>
  44 #include <SMESH_Gen.hxx>
  45 #include <SMESH_Group.hxx>
  46 #include <SMESH_Mesh.hxx>
  47 #include <SMESH_MeshEditor.hxx>
  48 #include <SMESH_MesherHelper.hxx>
  49 #include <StdMeshers_FaceSide.hxx>
  50 #include <StdMeshers_ViscousLayers2D.hxx>
  51
  52 #include <utilities.h>
  53
  54 #include <limits>
  55 #include <list>
  56 #include <vector>
  57 #include <set>
  58 #include <cstdlib>
  59
  60 // OPENCASCADE includes
  61 #include <BRep_Tool.hxx>
  62 #include <TopExp.hxx>
  63 #include <TopExp_Explorer.hxx>
  64 #include <TopoDS.hxx>
  65 #include <NCollection_Map.hxx>
  66
  67 #include <Geom_Surface.hxx>
  68 #include <Handle_Geom_Surface.hxx>
  69 #include <Geom2d_Curve.hxx>
  70 #include <Handle_Geom2d_Curve.hxx>
  71 #include <Geom_Curve.hxx>
  72 #include <Handle_Geom_Curve.hxx>
  73 #include <Handle_AIS_InteractiveObject.hxx>
  74 #include <TopoDS_Vertex.hxx>
  75 #include <TopoDS_Edge.hxx>
  76 #include <TopoDS_Wire.hxx>
  77 #include <TopoDS_Face.hxx>
  78
  79 #include <gp_Pnt2d.hxx>
  80 #include <gp_Pnt.hxx>
  81 #include <TopTools_IndexedMapOfShape.hxx>
  82 #include <TopoDS_Shape.hxx>
  83 #include <BRep_Builder.hxx>
  84 #include <BRepBuilderAPI_MakeVertex.hxx>
  85 #include <BRepTools.hxx>
  86
  87 #include <TopTools_DataMapOfShapeInteger.hxx>
  88 #include <GProp_GProps.hxx>
  89 #include <BRepGProp.hxx>
  90
  91 #ifndef WNT
  92 #include <fenv.h>
  93 #endif
  94
  95 #include <Standard_ErrorHandler.hxx>
  96 #include <GeomAPI_ProjectPointOnCurve.hxx>
  97 #include <GeomAPI_ProjectPointOnSurf.hxx>
  98 #include <gp_XY.hxx>
  99 #include <gp_XYZ.hxx>
 100 #include <TopTools_MapOfShape.hxx>
 101
 102 /* ==================================
 103  * ===========  PYTHON ==============
 104  * ==================================*/
 105
 106 typedef struct {
 107   PyObject_HEAD
 108   int softspace;
 109   std::string *out;
 110   } PyStdOut;
 111
 112 static void
 113 PyStdOut_dealloc(PyStdOut *self)
 114 {
 115   PyObject_Del(self);
 116 }
 117
 118 static PyObject *
 119 PyStdOut_write(PyStdOut *self, PyObject *args)
 120 {
 121   char *c;
 122   int l;
 123   if (!PyArg_ParseTuple(args, "t#:write",&c, &l))
 124     return NULL;
 125
 126   //std::cerr << c ;
 127   *(self->out)=*(self->out)+c;
 128
 129   Py_INCREF(Py_None);
 130   return Py_None;
 131 }
 132
 133 static PyMethodDef PyStdOut_methods[] = {
 134   {"write",  (PyCFunction)PyStdOut_write,  METH_VARARGS,
 135     PyDoc_STR("write(string) -> None")},
 136   {NULL,    NULL}   /* sentinel */
 137 };
 138
 139 static PyMemberDef PyStdOut_memberlist[] = {
 140   {(char*)"softspace", T_INT,  offsetof(PyStdOut, softspace), 0,
 141    (char*)"flag indicating that a space needs to be printed; used by print"},
 142   {NULL} /* Sentinel */
 143 };
 144
 145 static PyTypeObject PyStdOut_Type = {
 146   /* The ob_type field must be initialized in the module init function
 147    * to be portable to Windows without using C++. */
 148   PyObject_HEAD_INIT(NULL)
 149   0,                            /*ob_size*/
 150   "PyOut",                      /*tp_name*/
 151   sizeof(PyStdOut),             /*tp_basicsize*/
 152   0,                            /*tp_itemsize*/
 153   /* methods */
 154   (destructor)PyStdOut_dealloc, /*tp_dealloc*/
 155   0,                            /*tp_print*/
 156   0,                            /*tp_getattr*/
 157   0,                            /*tp_setattr*/
 158   0,                            /*tp_compare*/
 159   0,                            /*tp_repr*/
 160   0,                            /*tp_as_number*/
 161   0,                            /*tp_as_sequence*/
 162   0,                            /*tp_as_mapping*/
 163   0,                            /*tp_hash*/
 164   0,                            /*tp_call*/
 165   0,                            /*tp_str*/
 166   PyObject_GenericGetAttr,      /*tp_getattro*/
 167   /* softspace is writable:  we must supply tp_setattro */
 168   PyObject_GenericSetAttr,      /* tp_setattro */
 169   0,                            /*tp_as_buffer*/
 170   Py_TPFLAGS_DEFAULT,           /*tp_flags*/
 171   0,                            /*tp_doc*/
 172   0,                            /*tp_traverse*/
 173   0,                            /*tp_clear*/
 174   0,                            /*tp_richcompare*/
 175   0,                            /*tp_weaklistoffset*/
 176   0,                            /*tp_iter*/
 177   0,                            /*tp_iternext*/
 178   PyStdOut_methods,             /*tp_methods*/
 179   PyStdOut_memberlist,          /*tp_members*/
 180   0,                            /*tp_getset*/
 181   0,                            /*tp_base*/
 182   0,                            /*tp_dict*/
 183   0,                            /*tp_descr_get*/
 184   0,                            /*tp_descr_set*/
 185   0,                            /*tp_dictoffset*/
 186   0,                            /*tp_init*/
 187   0,                            /*tp_alloc*/
 188   0,                            /*tp_new*/
 189   0,                            /*tp_free*/
 190   0,                            /*tp_is_gc*/
 191 };
 192
 193 PyObject * newPyStdOut( std::string& out )
 194 {
 195   PyStdOut *self;
 196   self = PyObject_New(PyStdOut, &PyStdOut_Type);
 197   if (self == NULL)
 198     return NULL;
 199   self->softspace = 0;
 200   self->out=&out;
 201   return (PyObject*)self;
 202 }
 203
 204
 205 ////////////////////////END PYTHON///////////////////////////
 206
 207 //////////////////MY MAPS////////////////////////////////////////
 208 TopTools_IndexedMapOfShape FacesWithSizeMap;
 209 std::map<int,string> FaceId2SizeMap;
 210 TopTools_IndexedMapOfShape EdgesWithSizeMap;
 211 std::map<int,string> EdgeId2SizeMap;
 212 TopTools_IndexedMapOfShape VerticesWithSizeMap;
 213 std::map<int,string> VertexId2SizeMap;
 214
 215 std::map<int,PyObject*> FaceId2PythonSmp;
 216 std::map<int,PyObject*> EdgeId2PythonSmp;
 217 std::map<int,PyObject*> VertexId2PythonSmp;
 218
 219 std::map<int,std::vector<double> > FaceId2AttractorCoords;
 220 std::map<int,BLSURFPlugin_Attractor*> FaceId2ClassAttractor;
 221 std::map<int,BLSURFPlugin_Attractor*> FaceIndex2ClassAttractor;
 222
 223 TopTools_IndexedMapOfShape FacesWithEnforcedVertices;
 224 std::map< int, BLSURFPlugin_Hypothesis::TEnfVertexCoordsList > FaceId2EnforcedVertexCoords;
 225 std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexCoords > EnfVertexCoords2ProjVertex;
 226 std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList > EnfVertexCoords2EnfVertexList;
 227
 228 bool HasSizeMapOnFace=false;
 229 bool HasSizeMapOnEdge=false;
 230 bool HasSizeMapOnVertex=false;
 231 //bool HasAttractorOnFace=false;
 232
 233 //=============================================================================
 234 /*!
 235  *
 236  */
 237 //=============================================================================
 238
 239 BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF(int hypId, int studyId,
 240                                                SMESH_Gen* gen)
 241   : SMESH_2D_Algo(hypId, studyId, gen)
 242 {
 243   MESSAGE("BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF");
 244
 245   _name = "BLSURF";
 246   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
 247   _compatibleHypothesis.push_back(BLSURFPlugin_Hypothesis::GetHypType());
 248   _compatibleHypothesis.push_back(StdMeshers_ViscousLayers2D::GetHypType());
 249   _requireDiscreteBoundary = false;
 250   _onlyUnaryInput = false;
 251   _hypothesis = NULL;
 252   _supportSubmeshes = true;
 253
 254   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 255   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 256   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 257
 258   MESSAGE("studyid = " << _studyId);
 259
 260   myStudy = NULL;
 261   myStudy = aStudyMgr->GetStudyByID(_studyId);
 262   if (myStudy)
 263     MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
 264
 265   /* Initialize the Python interpreter */
 266   assert(Py_IsInitialized());
 267   PyGILState_STATE gstate;
 268   gstate = PyGILState_Ensure();
 269
 270   main_mod = NULL;
 271   main_mod = PyImport_AddModule("__main__");
 272
 273   main_dict = NULL;
 274   main_dict = PyModule_GetDict(main_mod);
 275
 276   PyRun_SimpleString("from math import *");
 277   PyGILState_Release(gstate);
 278
 279   FacesWithSizeMap.Clear();
 280   FaceId2SizeMap.clear();
 281   EdgesWithSizeMap.Clear();
 282   EdgeId2SizeMap.clear();
 283   VerticesWithSizeMap.Clear();
 284   VertexId2SizeMap.clear();
 285   FaceId2PythonSmp.clear();
 286   EdgeId2PythonSmp.clear();
 287   VertexId2PythonSmp.clear();
 288   FaceId2AttractorCoords.clear();
 289   FaceId2ClassAttractor.clear();
 290   FaceIndex2ClassAttractor.clear();
 291   FacesWithEnforcedVertices.Clear();
 292   FaceId2EnforcedVertexCoords.clear();
 293   EnfVertexCoords2ProjVertex.clear();
 294   EnfVertexCoords2EnfVertexList.clear();
 295
 296 #ifdef WITH_SMESH_CANCEL_COMPUTE
 297   _compute_canceled = false;
 298 #endif
 299 }
 300
 301 //=============================================================================
 302 /*!
 303  *
 304  */
 305 //=============================================================================
 306
 307 BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF()
 308 {
 309   MESSAGE("BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF");
 310 }
 311
 312
 313 //=============================================================================
 314 /*!
 315  *
 316  */
 317 //=============================================================================
 318
 319 bool BLSURFPlugin_BLSURF::CheckHypothesis
 320                          (SMESH_Mesh&                          aMesh,
 321                           const TopoDS_Shape&                  aShape,
 322                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 323 {
 324   _hypothesis        = NULL;
 325   _haveViscousLayers = false;
 326
 327   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 328   const SMESHDS_Hypothesis* theHyp;
 329
 330   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape,
 331                                                                   /*ignoreAuxiliary=*/false);
 332   aStatus = SMESH_Hypothesis::HYP_OK;
 333   if ( hyps.empty() )
 334   {
 335     return true;  // can work with no hypothesis
 336   }
 337
 338   for ( itl = hyps.begin(); itl != hyps.end(); ++itl )
 339   {
 340     theHyp = *itl;
 341     string hypName = theHyp->GetName();
 342     if ( hypName == BLSURFPlugin_Hypothesis::GetHypType() )
 343     {
 344       _hypothesis = static_cast<const BLSURFPlugin_Hypothesis*> (theHyp);
 345       ASSERT(_hypothesis);
 346       if ( _hypothesis->GetPhysicalMesh() == BLSURFPlugin_Hypothesis::DefaultSize &&
 347            _hypothesis->GetGeometricMesh() == BLSURFPlugin_Hypothesis::DefaultGeom )
 348         //  hphy_flag = 0 and hgeo_flag = 0 is not allowed (spec)
 349         aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
 350     }
 351     else if ( hypName == StdMeshers_ViscousLayers2D::GetHypType() )
 352     {
 353       _haveViscousLayers = true;
 354     }
 355     else
 356     {
 357       aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 358     }
 359   }
 360   return aStatus == SMESH_Hypothesis::HYP_OK;
 361 }
 362
 363 //=============================================================================
 364 /*!
 365  * Pass parameters to BLSURF
 366  */
 367 //=============================================================================
 368
 369 inline std::string to_string(double d)
 370 {
 371    std::ostringstream o;
 372    o << d;
 373    return o.str();
 374 }
 375
 376 inline std::string to_string_rel(double d)
 377 {
 378    std::ostringstream o;
 379    o << d;
 380    o << 'r';
 381    return o.str();
 382 }
 383
 384 inline std::string to_string(int i)
 385 {
 386    std::ostringstream o;
 387    o << i;
 388    return o.str();
 389 }
 390
 391 inline std::string to_string_rel(int i)
 392 {
 393    std::ostringstream o;
 394    o << i;
 395    o << 'r';
 396    return o.str();
 397 }
 398
 399 double _smp_phy_size;
 400 // #if BLSURF_VERSION_LONG >= "3.1.1"
 401 // //   sizemap_t *geo_sizemap_e, *geo_sizemap_f;
 402 //   sizemap_t *iso_sizemap_p, *iso_sizemap_e, *iso_sizemap_f;
 403 // //   sizemap_t *clean_geo_sizemap_e, *clean_geo_sizemap_f;
 404 //   sizemap_t *clean_iso_sizemap_p, *clean_iso_sizemap_e, *clean_iso_sizemap_f;
 405 // #endif
 406 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data);
 407 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data);
 408 status_t size_on_vertex(integer vertex_id, real *size, void *user_data);
 409
 410 typedef struct {
 411         gp_XY uv;
 412         gp_XYZ xyz;
 413 } projectionPoint;
 414 /////////////////////////////////////////////////////////
 415 projectionPoint getProjectionPoint(const TopoDS_Face& face, const gp_Pnt& point)
 416 {
 417   projectionPoint myPoint;
 418   Handle(Geom_Surface) surface = BRep_Tool::Surface(face);
 419   GeomAPI_ProjectPointOnSurf projector( point, surface );
 420   if ( !projector.IsDone() || projector.NbPoints()==0 )
 421     throw "getProjectionPoint: Can't project";
 422
 423   Quantity_Parameter u,v;
 424   projector.LowerDistanceParameters(u,v);
 425   myPoint.uv = gp_XY(u,v);
 426   gp_Pnt aPnt = projector.NearestPoint();
 427   myPoint.xyz = gp_XYZ(aPnt.X(),aPnt.Y(),aPnt.Z());
 428   //return gp_XY(u,v);
 429   return myPoint;
 430 }
 431 /////////////////////////////////////////////////////////
 432
 433 /////////////////////////////////////////////////////////
 434 double getT(const TopoDS_Edge& edge, const gp_Pnt& point)
 435 {
 436   Standard_Real f,l;
 437   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, f,l);
 438   GeomAPI_ProjectPointOnCurve projector( point, curve);
 439   if ( projector.NbPoints() == 0 )
 440     throw;
 441   return projector.LowerDistanceParameter();
 442 }
 443
 444 /////////////////////////////////////////////////////////
 445 TopoDS_Shape BLSURFPlugin_BLSURF::entryToShape(std::string entry)
 446 {
 447   MESSAGE("BLSURFPlugin_BLSURF::entryToShape "<<entry );
 448   GEOM::GEOM_Object_var aGeomObj;
 449   TopoDS_Shape S = TopoDS_Shape();
 450   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 451   SALOMEDS::GenericAttribute_var anAttr;
 452
 453   if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
 454     SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 455     CORBA::String_var aVal = anIOR->Value();
 456     CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
 457     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 458   }
 459   if ( !aGeomObj->_is_nil() )
 460     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 461   return S;
 462 }
 463
 464 void _createEnforcedVertexOnFace(TopoDS_Face faceShape, gp_Pnt aPnt, BLSURFPlugin_Hypothesis::TEnfVertex *enfVertex)
 465 {
 466   BLSURFPlugin_Hypothesis::TEnfVertexCoords enf_coords, coords, s_coords;
 467   enf_coords.clear();
 468   coords.clear();
 469   s_coords.clear();
 470
 471   // Get the (u,v) values of the enforced vertex on the face
 472   projectionPoint myPoint = getProjectionPoint(faceShape,aPnt);
 473
 474   MESSAGE("Enforced Vertex: " << aPnt.X() << ", " << aPnt.Y() << ", " << aPnt.Z());
 475   MESSAGE("Projected Vertex: " << myPoint.xyz.X() << ", " << myPoint.xyz.Y() << ", " << myPoint.xyz.Z());
 476   MESSAGE("Parametric coordinates: " << myPoint.uv.X() << ", " << myPoint.uv.Y() );
 477
 478   enf_coords.push_back(aPnt.X());
 479   enf_coords.push_back(aPnt.Y());
 480   enf_coords.push_back(aPnt.Z());
 481
 482   coords.push_back(myPoint.uv.X());
 483   coords.push_back(myPoint.uv.Y());
 484   coords.push_back(myPoint.xyz.X());
 485   coords.push_back(myPoint.xyz.Y());
 486   coords.push_back(myPoint.xyz.Z());
 487
 488   s_coords.push_back(myPoint.xyz.X());
 489   s_coords.push_back(myPoint.xyz.Y());
 490   s_coords.push_back(myPoint.xyz.Z());
 491
 492   // Save pair projected vertex / enf vertex
 493   MESSAGE("Storing pair projected vertex / enf vertex:");
 494   MESSAGE("("<< myPoint.xyz.X() << ", " << myPoint.xyz.Y() << ", " << myPoint.xyz.Z() <<") / (" << aPnt.X() << ", " << aPnt.Y() << ", " << aPnt.Z()<<")");
 495   EnfVertexCoords2ProjVertex[s_coords] = enf_coords;
 496   MESSAGE("Group name is: \"" << enfVertex->grpName << "\"");
 497   pair<BLSURFPlugin_Hypothesis::TEnfVertexList::iterator,bool> ret;
 498   BLSURFPlugin_Hypothesis::TEnfVertexList::iterator it;
 499   ret = EnfVertexCoords2EnfVertexList[s_coords].insert(enfVertex);
 500   if (ret.second == false) {
 501     it = ret.first;
 502     (*it)->grpName = enfVertex->grpName;
 503   }
 504
 505   int key = 0;
 506   if (! FacesWithEnforcedVertices.Contains(faceShape)) {
 507     key = FacesWithEnforcedVertices.Add(faceShape);
 508   }
 509   else {
 510     key = FacesWithEnforcedVertices.FindIndex(faceShape);
 511   }
 512
 513   // If a node is already created by an attractor, do not create enforced vertex
 514   int attractorKey = FacesWithSizeMap.FindIndex(faceShape);
 515   bool sameAttractor = false;
 516   if (attractorKey >= 0)
 517     if (FaceId2AttractorCoords.count(attractorKey) > 0)
 518       if (FaceId2AttractorCoords[attractorKey] == coords)
 519         sameAttractor = true;
 520
 521   if (FaceId2EnforcedVertexCoords.find(key) != FaceId2EnforcedVertexCoords.end()) {
 522     MESSAGE("Map of enf. vertex has key " << key)
 523     MESSAGE("Enf. vertex list size is: " << FaceId2EnforcedVertexCoords[key].size())
 524     if (! sameAttractor)
 525       FaceId2EnforcedVertexCoords[key].insert(coords); // there should be no redondant coords here (see std::set management)
 526     else
 527       MESSAGE("An attractor node is already defined: I don't add the enforced vertex");
 528     MESSAGE("New Enf. vertex list size is: " << FaceId2EnforcedVertexCoords[key].size())
 529   }
 530   else {
 531     MESSAGE("Map of enf. vertex has not key " << key << ": creating it")
 532     if (! sameAttractor) {
 533       BLSURFPlugin_Hypothesis::TEnfVertexCoordsList ens;
 534       ens.insert(coords);
 535       FaceId2EnforcedVertexCoords[key] = ens;
 536     }
 537     else
 538       MESSAGE("An attractor node is already defined: I don't add the enforced vertex");
 539   }
 540 }
 541
 542 /////////////////////////////////////////////////////////
 543 void BLSURFPlugin_BLSURF::createEnforcedVertexOnFace(TopoDS_Shape faceShape, BLSURFPlugin_Hypothesis::TEnfVertexList enfVertexList)
 544 {
 545   BLSURFPlugin_Hypothesis::TEnfVertex* enfVertex;
 546   gp_Pnt aPnt;
 547
 548   BLSURFPlugin_Hypothesis::TEnfVertexList::const_iterator enfVertexListIt = enfVertexList.begin();
 549
 550   for( ; enfVertexListIt != enfVertexList.end() ; ++enfVertexListIt ) {
 551     enfVertex = *enfVertexListIt;
 552     // Case of manual coords
 553     if (enfVertex->coords.size() != 0) {
 554       aPnt.SetCoord(enfVertex->coords[0],enfVertex->coords[1],enfVertex->coords[2]);
 555       _createEnforcedVertexOnFace( TopoDS::Face(faceShape),  aPnt, enfVertex);
 556     }
 557
 558     // Case of geom vertex coords
 559     if (enfVertex->geomEntry != "") {
 560       TopoDS_Shape GeomShape = entryToShape(enfVertex->geomEntry);
 561       TopAbs_ShapeEnum GeomType  = GeomShape.ShapeType();
 562        if (GeomType == TopAbs_VERTEX){
 563          aPnt = BRep_Tool::Pnt(TopoDS::Vertex(GeomShape));
 564          _createEnforcedVertexOnFace( TopoDS::Face(faceShape),  aPnt, enfVertex);
 565        }
 566        // Group Management
 567        if (GeomType == TopAbs_COMPOUND){
 568          for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 569            if (it.Value().ShapeType() == TopAbs_VERTEX){
 570              aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value()));
 571              _createEnforcedVertexOnFace( TopoDS::Face(faceShape),  aPnt, enfVertex);
 572            }
 573          }
 574        }
 575     }
 576   }
 577 }
 578
 579 /////////////////////////////////////////////////////////
 580 void createAttractorOnFace(TopoDS_Shape GeomShape, std::string AttractorFunction, double defaultSize)
 581 {
 582   MESSAGE("Attractor function: "<< AttractorFunction);
 583   double xa, ya, za; // Coordinates of attractor point
 584   double a, b;       // Attractor parameter
 585   double d = 0.;
 586   bool createNode=false; // To create a node on attractor projection
 587   int pos1, pos2;
 588   const char *sep = ";";
 589   // atIt->second has the following pattern:
 590   // ATTRACTOR(xa;ya;za;a;b;True|False;d)
 591   // where:
 592   // xa;ya;za : coordinates of  attractor
 593   // a        : desired size on attractor
 594   // b        : distance of influence of attractor
 595   // d        : distance until which the size remains constant
 596   //
 597   // We search the parameters in the string
 598   // xa
 599   pos1 = AttractorFunction.find(sep);
 600   if (pos1!=string::npos)
 601   xa = atof(AttractorFunction.substr(10, pos1-10).c_str());
 602   // ya
 603   pos2 = AttractorFunction.find(sep, pos1+1);
 604   if (pos2!=string::npos) {
 605   ya = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 606   pos1 = pos2;
 607   }
 608   // za
 609   pos2 = AttractorFunction.find(sep, pos1+1);
 610   if (pos2!=string::npos) {
 611   za = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 612   pos1 = pos2;
 613   }
 614   // a
 615   pos2 = AttractorFunction.find(sep, pos1+1);
 616   if (pos2!=string::npos) {
 617   a = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 618   pos1 = pos2;
 619   }
 620   // b
 621   pos2 = AttractorFunction.find(sep, pos1+1);
 622   if (pos2!=string::npos) {
 623   b = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 624   pos1 = pos2;
 625   }
 626   // createNode
 627   pos2 = AttractorFunction.find(sep, pos1+1);
 628   if (pos2!=string::npos) {
 629     string createNodeStr = AttractorFunction.substr(pos1+1, pos2-pos1-1);
 630     MESSAGE("createNode: " << createNodeStr);
 631     createNode = (AttractorFunction.substr(pos1+1, pos2-pos1-1) == "True");
 632     pos1=pos2;
 633   }
 634   // d
 635   pos2 = AttractorFunction.find(")");
 636   if (pos2!=string::npos) {
 637   d = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 638   }
 639
 640   // Get the (u,v) values of the attractor on the face
 641   projectionPoint myPoint = getProjectionPoint(TopoDS::Face(GeomShape),gp_Pnt(xa,ya,za));
 642   gp_XY uvPoint = myPoint.uv;
 643   gp_XYZ xyzPoint = myPoint.xyz;
 644   Standard_Real u0 = uvPoint.X();
 645   Standard_Real v0 = uvPoint.Y();
 646   Standard_Real x0 = xyzPoint.X();
 647   Standard_Real y0 = xyzPoint.Y();
 648   Standard_Real z0 = xyzPoint.Z();
 649   std::vector<double> coords;
 650   coords.push_back(u0);
 651   coords.push_back(v0);
 652   coords.push_back(x0);
 653   coords.push_back(y0);
 654   coords.push_back(z0);
 655   // We construct the python function
 656   ostringstream attractorFunctionStream;
 657   attractorFunctionStream << "def f(u,v): return ";
 658   attractorFunctionStream << defaultSize << "-(" << defaultSize <<"-" << a << ")";
 659   //attractorFunctionStream << "*exp(-((u-("<<u0<<"))*(u-("<<u0<<"))+(v-("<<v0<<"))*(v-("<<v0<<")))/(" << b << "*" << b <<"))";
 660   // rnc: make possible to keep the size constant until
 661   // a defined distance. Distance is expressed as the positiv part
 662   // of r-d where r is the distance to (u0,v0)
 663   attractorFunctionStream << "*exp(-(0.5*(sqrt((u-"<<u0<<")**2+(v-"<<v0<<")**2)-"<<d<<"+abs(sqrt((u-"<<u0<<")**2+(v-"<<v0<<")**2)-"<<d<<"))/(" << b << "))**2)";
 664
 665   MESSAGE("Python function for attractor:" << std::endl << attractorFunctionStream.str());
 666
 667   int key;
 668   if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 669     key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 670   }
 671   else {
 672     key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 673   }
 674   FaceId2SizeMap[key] =attractorFunctionStream.str();
 675   if (createNode) {
 676     MESSAGE("Creating node on ("<<x0<<","<<y0<<","<<z0<<")");
 677     FaceId2AttractorCoords[key] = coords;
 678   }
 679 //   // Test for new attractors
 680 //   gp_Pnt myP(xyzPoint);
 681 //   TopoDS_Vertex myV = BRepBuilderAPI_MakeVertex(myP);
 682 //   BLSURFPlugin_Attractor myAttractor(TopoDS::Face(GeomShape),myV,200);
 683 //   myAttractor.SetParameters(a, defaultSize, b, d);
 684 //   myAttractor.SetType(1);
 685 //   FaceId2ClassAttractor[key] = myAttractor;
 686 //   if(FaceId2ClassAttractor[key].GetFace().IsNull()){
 687 //     MESSAGE("face nulle ");
 688 //   }
 689 //   else
 690 //     MESSAGE("face OK");
 691 //
 692 //   if (FaceId2ClassAttractor[key].GetAttractorShape().IsNull()){
 693 //     MESSAGE("pas de point");
 694 //   }
 695 //   else
 696 //     MESSAGE("point OK");
 697 }
 698
 699 /////////////////////////////////////////////////////////
 700
 701 void BLSURFPlugin_BLSURF::SetParameters(
 702 // #if BLSURF_VERSION_LONG >= "3.1.1"
 703 //                                         cad_t *                          c,
 704 // #endif
 705                                         const BLSURFPlugin_Hypothesis* hyp,
 706                                         cadsurf_session_t *            css,
 707                                         precad_session_t *             pcs,
 708                                         SMESH_Mesh&                   mesh,
 709                                         bool *                  use_precad
 710                                        )
 711 {
 712   // rnc : Bug 1457
 713   // Clear map so that it is not stored in the algorithm with old enforced vertices in it
 714   EnfVertexCoords2EnfVertexList.clear();
 715
 716    double diagonal               = mesh.GetShapeDiagonalSize();
 717    double bbSegmentation         = _gen->GetBoundaryBoxSegmentation();
 718    int    _physicalMesh          = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
 719    int    _geometricMesh         = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
 720    double _phySize               = BLSURFPlugin_Hypothesis::GetDefaultPhySize(diagonal, bbSegmentation);
 721    bool   _phySizeRel            = BLSURFPlugin_Hypothesis::GetDefaultPhySizeRel();
 722    double _minSize               = BLSURFPlugin_Hypothesis::GetDefaultMinSize(diagonal);
 723    bool   _minSizeRel            = BLSURFPlugin_Hypothesis::GetDefaultMinSizeRel();
 724    double _maxSize               = BLSURFPlugin_Hypothesis::GetDefaultMaxSize(diagonal);
 725    bool   _maxSizeRel            = BLSURFPlugin_Hypothesis::GetDefaultMaxSizeRel();
 726    double _gradation             = BLSURFPlugin_Hypothesis::GetDefaultGradation();
 727    bool   _quadAllowed           = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
 728    double _angleMesh             = BLSURFPlugin_Hypothesis::GetDefaultAngleMesh();
 729    double _chordalError          = BLSURFPlugin_Hypothesis::GetDefaultChordalError(diagonal);
 730    bool   _anisotropic           = BLSURFPlugin_Hypothesis::GetDefaultAnisotropic();
 731    double _anisotropicRatio      = BLSURFPlugin_Hypothesis::GetDefaultAnisotropicRatio();
 732    bool   _removeTinyEdges       = BLSURFPlugin_Hypothesis::GetDefaultRemoveTinyEdges();
 733    double _tinyEdgeLength        = BLSURFPlugin_Hypothesis::GetDefaultTinyEdgeLength(diagonal);
 734    bool   _badElementRemoval     = BLSURFPlugin_Hypothesis::GetDefaultBadElementRemoval();
 735    double _badElementAspectRatio = BLSURFPlugin_Hypothesis::GetDefaultBadElementAspectRatio();
 736    bool   _optimizeMesh          = BLSURFPlugin_Hypothesis::GetDefaultOptimizeMesh();
 737    bool   _quadraticMesh         = BLSURFPlugin_Hypothesis::GetDefaultQuadraticMesh();
 738    int    _verb                  = BLSURFPlugin_Hypothesis::GetDefaultVerbosity();
 739    int    _topology              = BLSURFPlugin_Hypothesis::GetDefaultTopology();
 740
 741   // PreCAD
 742    int _precadMergeEdges         = BLSURFPlugin_Hypothesis::GetDefaultPreCADMergeEdges();
 743    int _precadProcess3DTopology  = BLSURFPlugin_Hypothesis::GetDefaultPreCADProcess3DTopology();
 744    int _precadDiscardInput       = BLSURFPlugin_Hypothesis::GetDefaultPreCADDiscardInput();
 745
 746
 747   if (hyp) {
 748     MESSAGE("BLSURFPlugin_BLSURF::SetParameters");
 749     _physicalMesh  = (int) hyp->GetPhysicalMesh();
 750     _geometricMesh = (int) hyp->GetGeometricMesh();
 751      if (hyp->GetPhySize() > 0)
 752        _phySize       = hyp->GetPhySize();
 753      _phySizeRel    = hyp->IsPhySizeRel();
 754      if (hyp->GetMinSize() > 0)
 755        _minSize       = hyp->GetMinSize();
 756      _minSizeRel    = hyp->IsMinSizeRel();
 757      if (hyp->GetMaxSize() > 0)
 758        _maxSize       = hyp->GetMaxSize();
 759      _maxSizeRel    = hyp->IsMaxSizeRel();
 760      if (hyp->GetGradation() > 0)
 761        _gradation     = hyp->GetGradation();
 762     _quadAllowed   = hyp->GetQuadAllowed();
 763      if (hyp->GetAngleMesh() > 0)
 764      _angleMesh     = hyp->GetAngleMesh();
 765      if (hyp->GetChordalError() > 0)
 766        _chordalError           = hyp->GetChordalError();
 767      _anisotropic            = hyp->GetAnisotropic();
 768      if (hyp->GetAnisotropicRatio() >= 0)
 769        _anisotropicRatio       = hyp->GetAnisotropicRatio();
 770      _removeTinyEdges        = hyp->GetRemoveTinyEdges();
 771      if (hyp->GetTinyEdgeLength() > 0)
 772        _tinyEdgeLength         = hyp->GetTinyEdgeLength();
 773      _badElementRemoval      = hyp->GetBadElementRemoval();
 774      if (hyp->GetBadElementAspectRatio() >= 0)
 775        _badElementAspectRatio  = hyp->GetBadElementAspectRatio();
 776      _optimizeMesh  = hyp->GetOptimizeMesh();
 777      _quadraticMesh = hyp->GetQuadraticMesh();
 778     _verb          = hyp->GetVerbosity();
 779      _topology      = (int) hyp->GetTopology();
 780      // PreCAD
 781      _precadMergeEdges = hyp->GetPreCADMergeEdges();
 782      _precadProcess3DTopology = hyp->GetPreCADProcess3DTopology();
 783      _precadDiscardInput = hyp->GetPreCADDiscardInput();
 784
 785     const BLSURFPlugin_Hypothesis::TOptionValues & opts = hyp->GetOptionValues();
 786     BLSURFPlugin_Hypothesis::TOptionValues::const_iterator opIt;
 787     for ( opIt = opts.begin(); opIt != opts.end(); ++opIt )
 788       if ( !opIt->second.empty() ) {
 789         MESSAGE("cadsurf_set_param(): " << opIt->first << " = " << opIt->second);
 790         cadsurf_set_param(css, opIt->first.c_str(), opIt->second.c_str());
 791       }
 792
 793     const BLSURFPlugin_Hypothesis::TOptionValues & preCADopts = hyp->GetPreCADOptionValues();
 794     for ( opIt = preCADopts.begin(); opIt != preCADopts.end(); ++opIt )
 795       if ( !opIt->second.empty() ) {
 796         if (_topology == BLSURFPlugin_Hypothesis::PreCAD) {
 797           MESSAGE("precad_set_param(): " << opIt->first << " = " << opIt->second);
 798           precad_set_param(pcs, opIt->first.c_str(), opIt->second.c_str());
 799         }
 800       }
 801   }
 802 //   else {
 803 //     //0020968: EDF1545 SMESH: Problem in the creation of a mesh group on geometry
 804 //     // GetDefaultPhySize() sometimes leads to computation failure
 805 //     // GDD 26/07/2012 From Distene documentation, global physical size default value = diag/100
 806 //     _phySize = BLSURFPlugin_Hypothesis::GetDefaultPhySize(diagonal);
 807 //     _minSize = BLSURFPlugin_Hypothesis::GetDefaultMinSize(diagonal);
 808 //     _maxSize = BLSURFPlugin_Hypothesis::GetDefaultMaxSize(diagonal);
 809 //     _chordalError = BLSURFPlugin_Hypothesis::GetDefaultChordalError(diagonal);
 810 //     _tinyEdgeLength = BLSURFPlugin_Hypothesis::GetDefaultTinyEdgeLength(diagonal);
 811 //     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
 812 //   }
 813
 814   // PreCAD
 815   if (_topology == BLSURFPlugin_Hypothesis::PreCAD) {
 816     *use_precad = true;
 817     precad_set_param(pcs, "verbose",                to_string(_verb).c_str());
 818     precad_set_param(pcs, "merge_edges",            _precadMergeEdges ? "1" : "0");
 819     precad_set_param(pcs, "process_3d_topology",    _precadProcess3DTopology ? "1" : "0");
 820     precad_set_param(pcs, "discard_input_topology", _precadDiscardInput ? "1" : "0");
 821   }
 822
 823    bool useGradation = false;
 824    switch (_physicalMesh)
 825    {
 826      case BLSURFPlugin_Hypothesis::PhysicalGlobalSize:
 827        cadsurf_set_param(css, "physical_size_mode", "global");
 828        cadsurf_set_param(css, "global_physical_size", _phySizeRel ? to_string_rel(_phySize).c_str() : to_string(_phySize).c_str());
 829        break;
 830      case BLSURFPlugin_Hypothesis::PhysicalLocalSize:
 831        cadsurf_set_param(css, "physical_size_mode", "local");
 832        cadsurf_set_param(css, "global_physical_size", _phySizeRel ? to_string_rel(_phySize).c_str() : to_string(_phySize).c_str());
 833        useGradation = true;
 834        break;
 835      default:
 836        cadsurf_set_param(css, "physical_size_mode", "none");
 837    }
 838
 839    switch (_geometricMesh)
 840    {
 841      case BLSURFPlugin_Hypothesis::GeometricalGlobalSize:
 842        cadsurf_set_param(css, "geometric_size_mode", "global");
 843        cadsurf_set_param(css, "geometric_approximation", to_string(_angleMesh).c_str());
 844        cadsurf_set_param(css, "chordal_error", to_string(_chordalError).c_str());
 845        useGradation = true;
 846        break;
 847      case BLSURFPlugin_Hypothesis::GeometricalLocalSize:
 848        cadsurf_set_param(css, "geometric_size_mode", "local");
 849        cadsurf_set_param(css, "geometric_approximation", to_string(_angleMesh).c_str());
 850        cadsurf_set_param(css, "chordal_error", to_string(_chordalError).c_str());
 851        useGradation = true;
 852        break;
 853      default:
 854        cadsurf_set_param(css, "geometric_size_mode", "none");
 855    }
 856
 857    cadsurf_set_param(css, "minsize",                           _minSizeRel ? to_string_rel(_minSize).c_str() : to_string(_minSize).c_str());
 858    cadsurf_set_param(css, "maxsize",                           _maxSizeRel ? to_string_rel(_maxSize).c_str() : to_string(_maxSize).c_str());
 859    if ( useGradation )
 860      cadsurf_set_param(css, "gradation",                         to_string(_gradation).c_str());
 861    cadsurf_set_param(css, "element_generation",                _quadAllowed ? "quad_dominant" : "triangle");
 862
 863
 864    cadsurf_set_param(css, "metric",                            _anisotropic ? "anisotropic" : "isotropic");
 865    if ( _anisotropic )
 866      cadsurf_set_param(css, "anisotropic_ratio",                 to_string(_anisotropicRatio).c_str());
 867    cadsurf_set_param(css, "remove_tiny_edges",                 _removeTinyEdges ? "1" : "0");
 868    if ( _removeTinyEdges )
 869      cadsurf_set_param(css, "tiny_edge_length",                  to_string(_tinyEdgeLength).c_str());
 870    cadsurf_set_param(css, "force_bad_surface_element_removal", _badElementRemoval ? "1" : "0");
 871    if ( _badElementRemoval )
 872      cadsurf_set_param(css, "bad_surface_element_aspect_ratio",  to_string(_badElementAspectRatio).c_str());
 873    cadsurf_set_param(css, "optimisation",                      _optimizeMesh ? "yes" : "no");
 874    cadsurf_set_param(css, "element_order",                     _quadraticMesh ? "quadratic" : "linear");
 875    cadsurf_set_param(css, "verbose",                           to_string(_verb).c_str());
 876
 877    _smp_phy_size = _phySizeRel ? _phySize*diagonal : _phySize;
 878    std::cout << "_smp_phy_size = " << _smp_phy_size << std::endl;
 879
 880    if (_physicalMesh == BLSURFPlugin_Hypothesis::PhysicalLocalSize){
 881     TopoDS_Shape GeomShape;
 882     TopoDS_Shape AttShape;
 883     TopAbs_ShapeEnum GeomType;
 884     //
 885     // Standard Size Maps
 886     //
 887     MESSAGE("Setting a Size Map");
 888     const BLSURFPlugin_Hypothesis::TSizeMap sizeMaps = BLSURFPlugin_Hypothesis::GetSizeMapEntries(hyp);
 889     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator smIt = sizeMaps.begin();
 890     for ( ; smIt != sizeMaps.end(); ++smIt ) {
 891       if ( !smIt->second.empty() ) {
 892         MESSAGE("cadsurf_set_sizeMap(): " << smIt->first << " = " << smIt->second);
 893         GeomShape = entryToShape(smIt->first);
 894         GeomType  = GeomShape.ShapeType();
 895         MESSAGE("Geomtype is " << GeomType);
 896         int key = -1;
 897         // Group Management
 898         if (GeomType == TopAbs_COMPOUND){
 899           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 900             // Group of faces
 901             if (it.Value().ShapeType() == TopAbs_FACE){
 902               HasSizeMapOnFace = true;
 903               if (! FacesWithSizeMap.Contains(TopoDS::Face(it.Value()))) {
 904                 key = FacesWithSizeMap.Add(TopoDS::Face(it.Value()));
 905               }
 906               else {
 907                 key = FacesWithSizeMap.FindIndex(TopoDS::Face(it.Value()));
 908 //                 MESSAGE("Face with key " << key << " already in map");
 909               }
 910               FaceId2SizeMap[key] = smIt->second;
 911             }
 912             // Group of edges
 913             if (it.Value().ShapeType() == TopAbs_EDGE){
 914               HasSizeMapOnEdge = true;
 915               HasSizeMapOnFace = true;
 916               if (! EdgesWithSizeMap.Contains(TopoDS::Edge(it.Value()))) {
 917                 key = EdgesWithSizeMap.Add(TopoDS::Edge(it.Value()));
 918               }
 919               else {
 920                 key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(it.Value()));
 921 //                 MESSAGE("Edge with key " << key << " already in map");
 922               }
 923               EdgeId2SizeMap[key] = smIt->second;
 924             }
 925             // Group of vertices
 926             if (it.Value().ShapeType() == TopAbs_VERTEX){
 927               HasSizeMapOnVertex = true;
 928               HasSizeMapOnEdge = true;
 929               HasSizeMapOnFace = true;
 930               if (! VerticesWithSizeMap.Contains(TopoDS::Vertex(it.Value()))) {
 931                 key = VerticesWithSizeMap.Add(TopoDS::Vertex(it.Value()));
 932               }
 933               else {
 934                 key = VerticesWithSizeMap.FindIndex(TopoDS::Vertex(it.Value()));
 935                 MESSAGE("Group of vertices with key " << key << " already in map");
 936               }
 937               MESSAGE("Group of vertices with key " << key << " has a size map: " << smIt->second);
 938               VertexId2SizeMap[key] = smIt->second;
 939             }
 940           }
 941         }
 942         // Single face
 943         if (GeomType == TopAbs_FACE){
 944           HasSizeMapOnFace = true;
 945           if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 946             key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 947           }
 948           else {
 949             key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 950 //             MESSAGE("Face with key " << key << " already in map");
 951           }
 952           FaceId2SizeMap[key] = smIt->second;
 953         }
 954         // Single edge
 955         if (GeomType == TopAbs_EDGE){
 956           HasSizeMapOnEdge = true;
 957           HasSizeMapOnFace = true;
 958           if (! EdgesWithSizeMap.Contains(TopoDS::Edge(GeomShape))) {
 959             key = EdgesWithSizeMap.Add(TopoDS::Edge(GeomShape));
 960           }
 961           else {
 962             key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(GeomShape));
 963 //             MESSAGE("Edge with key " << key << " already in map");
 964           }
 965           EdgeId2SizeMap[key] = smIt->second;
 966         }
 967         // Single vertex
 968         if (GeomType == TopAbs_VERTEX){
 969           HasSizeMapOnVertex = true;
 970           HasSizeMapOnEdge   = true;
 971           HasSizeMapOnFace   = true;
 972           if (! VerticesWithSizeMap.Contains(TopoDS::Vertex(GeomShape))) {
 973             key = VerticesWithSizeMap.Add(TopoDS::Vertex(GeomShape));
 974           }
 975           else {
 976             key = VerticesWithSizeMap.FindIndex(TopoDS::Vertex(GeomShape));
 977              MESSAGE("Vertex with key " << key << " already in map");
 978           }
 979           MESSAGE("Vertex with key " << key << " has a size map: " << smIt->second);
 980           VertexId2SizeMap[key] = smIt->second;
 981         }
 982       }
 983     }
 984
 985     //
 986     // Attractors
 987     //
 988     // TODO appeler le constructeur des attracteurs directement ici
 989     MESSAGE("Setting Attractors");
 990 //     if ( !_phySizeRel ) {
 991       const BLSURFPlugin_Hypothesis::TSizeMap attractors = BLSURFPlugin_Hypothesis::GetAttractorEntries(hyp);
 992       BLSURFPlugin_Hypothesis::TSizeMap::const_iterator atIt = attractors.begin();
 993       for ( ; atIt != attractors.end(); ++atIt ) {
 994         if ( !atIt->second.empty() ) {
 995           MESSAGE("cadsurf_set_attractor(): " << atIt->first << " = " << atIt->second);
 996           GeomShape = entryToShape(atIt->first);
 997           GeomType  = GeomShape.ShapeType();
 998           // Group Management
 999           if (GeomType == TopAbs_COMPOUND){
1000             for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
1001               if (it.Value().ShapeType() == TopAbs_FACE){
1002                 HasSizeMapOnFace = true;
1003                 createAttractorOnFace(it.Value(), atIt->second, _phySizeRel ? _phySize*diagonal : _phySize);
1004               }
1005             }
1006           }
1007
1008           if (GeomType == TopAbs_FACE){
1009             HasSizeMapOnFace = true;
1010             createAttractorOnFace(GeomShape, atIt->second, _phySizeRel ? _phySize*diagonal : _phySize);
1011           }
1012   /*
1013           if (GeomType == TopAbs_EDGE){
1014             HasSizeMapOnEdge = true;
1015             HasSizeMapOnFace = true;
1016           EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(IntegerLast())] = atIt->second;
1017           }
1018           if (GeomType == TopAbs_VERTEX){
1019             HasSizeMapOnVertex = true;
1020             HasSizeMapOnEdge   = true;
1021             HasSizeMapOnFace   = true;
1022           VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(IntegerLast())] = atIt->second;
1023           }
1024   */
1025         }
1026       }
1027 //     }
1028 //     else
1029 //       MESSAGE("Impossible to create the attractors when the physical size is relative");
1030
1031     // Class Attractors
1032     // temporary commented out for testing
1033     // TODO
1034     //  - Fill in the BLSURFPlugin_Hypothesis::TAttractorMap map in the hypothesis
1035     //  - Uncomment and complete this part to construct the attractors from the attractor shape and the passed parameters on each face of the map
1036     //  - To do this use the public methodss: SetParameters(several double parameters) and SetType(int type)
1037     //  OR, even better:
1038     //  - Construct the attractors with an empty dist. map in the hypothesis
1039     //  - build the map here for each face with an attractor set and only if the attractor shape as changed since the last call to _buildmap()
1040     //  -> define a bool _mapbuilt in the class that is set to false by default and set to true when calling _buildmap()  OK
1041
1042     const BLSURFPlugin_Hypothesis::TAttractorMap class_attractors = BLSURFPlugin_Hypothesis::GetClassAttractorEntries(hyp);
1043     int key=-1;
1044     BLSURFPlugin_Hypothesis::TAttractorMap::const_iterator AtIt = class_attractors.begin();
1045     for ( ; AtIt != class_attractors.end(); ++AtIt ) {
1046       if ( !AtIt->second->Empty() ) {
1047        // MESSAGE("cadsurf_set_attractor(): " << AtIt->first << " = " << AtIt->second);
1048         GeomShape = entryToShape(AtIt->first);
1049         AttShape = AtIt->second->GetAttractorShape();
1050         GeomType  = GeomShape.ShapeType();
1051         // Group Management
1052 //         if (GeomType == TopAbs_COMPOUND){
1053 //           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
1054 //             if (it.Value().ShapeType() == TopAbs_FACE){
1055 //               HasAttractorOnFace = true;
1056 //               myAttractor = BLSURFPluginAttractor(GeomShape, AttShape);
1057 //             }
1058 //           }
1059 //         }
1060
1061         if (GeomType == TopAbs_FACE
1062           && (AttShape.ShapeType() == TopAbs_VERTEX
1063            || AttShape.ShapeType() == TopAbs_EDGE
1064            || AttShape.ShapeType() == TopAbs_WIRE
1065            || AttShape.ShapeType() == TopAbs_COMPOUND) ){
1066             HasSizeMapOnFace = true;
1067
1068             if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape)) ) {
1069                 key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape) );
1070             }
1071             else {
1072               key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
1073 //                 MESSAGE("Face with key " << key << " already in map");
1074             }
1075
1076             FaceId2ClassAttractor[key] = AtIt->second;
1077         }
1078         else{
1079           MESSAGE("Wrong shape type !!")
1080         }
1081
1082       }
1083     }
1084
1085
1086     //
1087     // Enforced Vertices
1088     //
1089     MESSAGE("Setting Enforced Vertices");
1090     const BLSURFPlugin_Hypothesis::TFaceEntryEnfVertexListMap entryEnfVertexListMap = BLSURFPlugin_Hypothesis::GetAllEnforcedVerticesByFace(hyp);
1091     BLSURFPlugin_Hypothesis::TFaceEntryEnfVertexListMap::const_iterator enfIt = entryEnfVertexListMap.begin();
1092     for ( ; enfIt != entryEnfVertexListMap.end(); ++enfIt ) {
1093       if ( !enfIt->second.empty() ) {
1094         GeomShape = entryToShape(enfIt->first);
1095         GeomType  = GeomShape.ShapeType();
1096         // Group Management
1097         if (GeomType == TopAbs_COMPOUND){
1098           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
1099             if (it.Value().ShapeType() == TopAbs_FACE){
1100               HasSizeMapOnFace = true;
1101               createEnforcedVertexOnFace(it.Value(), enfIt->second);
1102             }
1103           }
1104         }
1105
1106         if (GeomType == TopAbs_FACE){
1107           HasSizeMapOnFace = true;
1108           createEnforcedVertexOnFace(GeomShape, enfIt->second);
1109         }
1110       }
1111     }
1112
1113     // Internal vertices
1114     bool useInternalVertexAllFaces = BLSURFPlugin_Hypothesis::GetInternalEnforcedVertexAllFaces(hyp);
1115     if (useInternalVertexAllFaces) {
1116       std::string grpName = BLSURFPlugin_Hypothesis::GetInternalEnforcedVertexAllFacesGroup(hyp);
1117       MESSAGE("Setting Internal Enforced Vertices");
1118       GeomShape = mesh.GetShapeToMesh();
1119       gp_Pnt aPnt;
1120       TopExp_Explorer exp (GeomShape, TopAbs_FACE);
1121       for (; exp.More(); exp.Next()){
1122         MESSAGE("Iterating shapes. Shape type is " << exp.Current().ShapeType());
1123         TopExp_Explorer exp_face (exp.Current(), TopAbs_VERTEX, TopAbs_EDGE);
1124         for (; exp_face.More(); exp_face.Next())
1125         {
1126           // Get coords of vertex
1127           // Check if current coords is already in enfVertexList
1128           // If coords not in enfVertexList, add new enfVertex
1129           aPnt = BRep_Tool::Pnt(TopoDS::Vertex(exp_face.Current()));
1130           MESSAGE("Found vertex on face at " << aPnt.X() <<", "<<aPnt.Y()<<", "<<aPnt.Z());
1131           BLSURFPlugin_Hypothesis::TEnfVertex* enfVertex = new BLSURFPlugin_Hypothesis::TEnfVertex();
1132           enfVertex->coords.push_back(aPnt.X());
1133           enfVertex->coords.push_back(aPnt.Y());
1134           enfVertex->coords.push_back(aPnt.Z());
1135           enfVertex->name = "";
1136           enfVertex->faceEntries.clear();
1137           enfVertex->geomEntry = "";
1138           enfVertex->grpName = grpName;
1139           enfVertex->vertex = TopoDS::Vertex( exp_face.Current() );
1140           _createEnforcedVertexOnFace( TopoDS::Face(exp.Current()),  aPnt, enfVertex);
1141           HasSizeMapOnFace = true;
1142         }
1143       }
1144     }
1145
1146     MESSAGE("Setting Size Map on FACES ");
1147 // #if BLSURF_VERSION_LONG < "3.1.1"
1148     cadsurf_data_set_sizemap_iso_cad_face(css, size_on_surface, &_smp_phy_size);
1149 // #else
1150 //     if (*use_precad)
1151 //       iso_sizemap_f = sizemap_new(c, distene_sizemap_type_iso_cad_face, (void *)size_on_surface, NULL);
1152 //     else
1153 //       clean_iso_sizemap_f = sizemap_new(c, distene_sizemap_type_iso_cad_face, (void *)size_on_surface, NULL);
1154 // #endif
1155
1156     if (HasSizeMapOnEdge){
1157       MESSAGE("Setting Size Map on EDGES ");
1158 // #if BLSURF_VERSION_LONG < "3.1.1"
1159       cadsurf_data_set_sizemap_iso_cad_edge(css, size_on_edge, &_smp_phy_size);
1160 // #else
1161 //       if (*use_precad)
1162 //         iso_sizemap_e = sizemap_new(c, distene_sizemap_type_iso_cad_edge, (void *)size_on_edge, NULL);
1163 //       else
1164 //         clean_iso_sizemap_e = sizemap_new(c, distene_sizemap_type_iso_cad_edge, (void *)size_on_edge, NULL);
1165 // #endif
1166     }
1167     if (HasSizeMapOnVertex){
1168       MESSAGE("Setting Size Map on VERTICES ");
1169 // #if BLSURF_VERSION_LONG < "3.1.1"
1170       cadsurf_data_set_sizemap_iso_cad_point(css, size_on_vertex, &_smp_phy_size);
1171 // #else
1172 //       if (*use_precad)
1173 //         iso_sizemap_p = sizemap_new(c, distene_sizemap_type_iso_cad_point, (void *)size_on_vertex, NULL);
1174 //       else
1175 //         clean_iso_sizemap_p = sizemap_new(c, distene_sizemap_type_iso_cad_point, (void *)size_on_vertex, NULL);
1176 // #endif
1177     }
1178   }
1179 }
1180
1181 namespace
1182 {
1183   // --------------------------------------------------------------------------
1184   /*!
1185    * \brief Class correctly terminating usage of BLSURF library at destruction
1186    */
1187   class BLSURF_Cleaner
1188   {
1189     context_t *       _ctx;
1190     cadsurf_session_t* _css;
1191     cad_t *           _cad;
1192     dcad_t *          _dcad;
1193   public:
1194     BLSURF_Cleaner(context_t *       ctx,
1195                    cadsurf_session_t* css,
1196                    cad_t *           cad,
1197                    dcad_t *          dcad)
1198       : _ctx ( ctx  ),
1199         _css ( css  ),
1200         _cad ( cad  ),
1201         _dcad( dcad )
1202     {
1203     }
1204     ~BLSURF_Cleaner()
1205     {
1206       Clean( /*exceptContext=*/false );
1207     }
1208     void Clean(const bool exceptContext)
1209     {
1210       if ( _css )
1211       {
1212         cadsurf_session_delete(_css); _css = 0;
1213
1214         // #if BLSURF_VERSION_LONG >= "3.1.1"
1215         // //     if(geo_sizemap_e)
1216         // //       distene_sizemap_delete(geo_sizemap_e);
1217         // //     if(geo_sizemap_f)
1218         // //       distene_sizemap_delete(geo_sizemap_f);
1219         //     if(iso_sizemap_p)
1220         //       distene_sizemap_delete(iso_sizemap_p);
1221         //     if(iso_sizemap_e)
1222         //       distene_sizemap_delete(iso_sizemap_e);
1223         //     if(iso_sizemap_f)
1224         //       distene_sizemap_delete(iso_sizemap_f);
1225         //
1226         // //     if(clean_geo_sizemap_e)
1227         // //       distene_sizemap_delete(clean_geo_sizemap_e);
1228         // //     if(clean_geo_sizemap_f)
1229         // //       distene_sizemap_delete(clean_geo_sizemap_f);
1230         //     if(clean_iso_sizemap_p)
1231         //       distene_sizemap_delete(clean_iso_sizemap_p);
1232         //     if(clean_iso_sizemap_e)
1233         //       distene_sizemap_delete(clean_iso_sizemap_e);
1234         //     if(clean_iso_sizemap_f)
1235         //       distene_sizemap_delete(clean_iso_sizemap_f);
1236         // #endif
1237
1238         cad_delete(_cad); _cad = 0;
1239         dcad_delete(_dcad); _dcad = 0;
1240         if ( !exceptContext )
1241         {
1242           context_delete(_ctx); _ctx = 0;
1243         }
1244       }
1245     }
1246   };
1247
1248   // --------------------------------------------------------------------------
1249   /*!
1250    * \brief A curve of internal boundary of viscous layer
1251    */
1252   class TProxyPCurve : public Geom2d_Curve
1253   {
1254     const UVPtStructVec& _noDataVec;
1255     Handle(Geom2d_Curve) _pcurve;
1256     double               _sign;
1257
1258     int locateU( const double U, double& r ) const
1259     {
1260       const double eps = 1e-10;
1261       double normU =
1262         ( U - _noDataVec[0].param ) / ( _noDataVec.back().param - _noDataVec[0].param );
1263       int i = int( normU * ( _noDataVec.size() - 1 ));
1264       while ( _noDataVec[ i ].param > U + eps && i != 0 )
1265         --i;
1266       while ( i+1 < _noDataVec.size() && _noDataVec[ i+1 ].param < U - eps )
1267         ++i;
1268       //cout << "U " << U << " normU " << normU << " i " << i << " _noDataVec.size() " << _noDataVec.size() << endl;
1269       r = 0;
1270       if ( normU > eps && normU < 1.- eps && i+1 < _noDataVec.size() )
1271         r = ( U - _noDataVec[i].param ) /
1272           ( _noDataVec[i+1].param - _noDataVec[i].param );
1273       return i;
1274     }
1275   public:
1276     TProxyPCurve( const SMESH_ProxyMesh::SubMesh* edgeSM,
1277                   Handle(Geom2d_Curve)            pcurve)
1278       : _noDataVec( edgeSM->GetUVPtStructVec() ), _pcurve( pcurve )
1279     {
1280       _sign = ( _noDataVec.back().param > _noDataVec[0].param ? +1 : -1 );
1281     }
1282     void D0(const Standard_Real U, gp_Pnt2d& P) const
1283     {
1284       double r;
1285       int i = locateU( U, r );
1286       P.SetCoord( _noDataVec[i].u, _noDataVec[i].v );
1287       if ( r > 0 )
1288         P = ( 1-r ) * P.XY() + r * gp_XY( _noDataVec[i+1].u, _noDataVec[i+1].v );
1289       // cout << "U " << U << " \ti,r = ( "<< i << ", "<< r << " )\t P ( " << P.X() << ", " << P.Y() <<" )" << endl;
1290     }
1291     gp_Vec2d DN(const Standard_Real U,const Standard_Integer N) const
1292     {
1293       double r;
1294       int i    = locateU( U, r );
1295       double u = _noDataVec[i].param;
1296       if ( r > 0 )
1297         u = ( 1-r ) * u + r * _noDataVec[i+1].param;
1298       gp_Vec2d dn = _pcurve->DN( u, N );
1299       //cout << "U " << U << " D" << N << " " << dn.X() << ", " << dn.Y() << end;
1300       return dn;
1301     }
1302     // virtual methods not used by BLSURF plug-in
1303     void Reverse()                                         {}
1304     Standard_Real ReversedParameter(const double U ) const { return 1 - U; }
1305     Standard_Real FirstParameter() const                   { return 0;}
1306     Standard_Real LastParameter() const                    { return 1;}
1307     Standard_Boolean IsClosed() const                      { return false; }
1308     Standard_Boolean IsPeriodic() const                    { return false; }
1309     GeomAbs_Shape Continuity() const                       { return GeomAbs_G2; }
1310     Standard_Boolean IsCN(const int N) const               { return false; }
1311     void Transform(const gp_Trsf2d&)                       {}
1312     Handle_Geom2d_Geometry Copy() const                    {return Handle_Geom2d_Geometry();}
1313     void D1(const Standard_Real U,gp_Pnt2d& P,gp_Vec2d& V1) const {}
1314     void D2(const Standard_Real U,gp_Pnt2d& P,gp_Vec2d& V1,gp_Vec2d& V2) const {}
1315     void D3(const Standard_Real U,gp_Pnt2d& P,gp_Vec2d& V1,gp_Vec2d& V2,gp_Vec2d& V3) const {}
1316
1317     //DEFINE_STANDARD_RTTI(TProxyPCurve);
1318   };
1319   DEFINE_STANDARD_HANDLE(TProxyPCurve,Geom2d_Curve);
1320   // --------------------------------------------------------------------------
1321   // comparator to sort nodes and sub-meshes
1322   struct ShapeTypeCompare
1323   {
1324     // sort nodes by position in the following order:
1325     // SMDS_TOP_FACE=2, SMDS_TOP_EDGE=1, SMDS_TOP_VERTEX=0, SMDS_TOP_3DSPACE=3
1326     int operator()( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 ) const
1327     {
1328       SMDS_TypeOfPosition pos1 = n1->GetPosition()->GetTypeOfPosition();
1329       SMDS_TypeOfPosition pos2 = n2->GetPosition()->GetTypeOfPosition();
1330       if ( pos1 == pos2 ) return 0;
1331       if ( pos1 < pos2 || pos1 == SMDS_TOP_3DSPACE ) return 1;
1332       return -1;
1333     }
1334     // sort sub-meshes in order: EDGE, VERTEX
1335     bool operator()( const SMESHDS_SubMesh* s1, const SMESHDS_SubMesh* s2 ) const
1336     {
1337       int isVertex1 = ( s1 && s1->NbElements() == 0 );
1338       int isVertex2 = ( s2 && s2->NbElements() == 0 );
1339       if ( isVertex1 == isVertex2 )
1340         return s1 < s2;
1341       return isVertex1 < isVertex2;
1342     }
1343   };
1344
1345   //================================================================================
1346   /*!
1347    * \brief Fills groups on nodes to be merged
1348    */
1349   //================================================================================
1350
1351   void getNodeGroupsToMerge( const SMESHDS_SubMesh*                smDS,
1352                              const TopoDS_Shape&                   shape,
1353                              SMESH_MeshEditor::TListOfListOfNodes& nodeGroupsToMerge)
1354   {
1355     SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
1356     switch ( shape.ShapeType() )
1357     {
1358     case TopAbs_VERTEX: {
1359       std::list< const SMDS_MeshNode* > nodes;
1360       while ( nIt->more() )
1361         nodes.push_back( nIt->next() );
1362       if ( nodes.size() > 1 )
1363         nodeGroupsToMerge.push_back( nodes );
1364       break;
1365     }
1366     case TopAbs_EDGE: {
1367       std::multimap< double, const SMDS_MeshNode* > u2node;
1368       const SMDS_EdgePosition* ePos;
1369       while ( nIt->more() )
1370       {
1371         const SMDS_MeshNode* n = nIt->next();
1372         if (( ePos = dynamic_cast< const SMDS_EdgePosition* >( n->GetPosition() )))
1373           u2node.insert( make_pair( ePos->GetUParameter(), n ));
1374       }
1375       if ( u2node.size() < 2 ) return;
1376
1377       double tol = (( u2node.rbegin()->first - u2node.begin()->first ) / 20.) / u2node.size();
1378       std::multimap< double, const SMDS_MeshNode* >::iterator un2, un1;
1379       for ( un2 = u2node.begin(), un1 = un2++; un2 != u2node.end(); un1 = un2++ )
1380       {
1381         if (( un2->first - un1->first ) <= tol )
1382         {
1383           std::list< const SMDS_MeshNode* > nodes;
1384           nodes.push_back( un1->second );
1385           while (( un2->first - un1->first ) <= tol )
1386           {
1387             nodes.push_back( un2->second );
1388             if ( ++un2 == u2node.end()) {
1389               --un2;
1390               break;
1391             }
1392           }
1393           // make nodes created on the boundary of viscous layer replace nodes created
1394           // by BLSURF as their SMDS_Position is more correct
1395           nodes.sort( ShapeTypeCompare() );
1396           nodeGroupsToMerge.push_back( nodes );
1397         }
1398       }
1399       break;
1400     }
1401     default: ;
1402     }
1403     // SMESH_MeshEditor::TListOfListOfNodes::const_iterator nll = nodeGroupsToMerge.begin();
1404     // for ( ; nll != nodeGroupsToMerge.end(); ++nll )
1405     // {
1406     //   cout << "Merge ";
1407     //   const std::list< const SMDS_MeshNode* >& nl = *nll;
1408     //   std::list< const SMDS_MeshNode* >::const_iterator nIt = nl.begin();
1409     //   for ( ; nIt != nl.end(); ++nIt )
1410     //     cout << (*nIt) << " ";
1411     //   cout << endl;
1412     // }
1413     // cout << endl;
1414   }
1415
1416 } // namespace
1417
1418 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
1419 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
1420                   real *duu, real *duv, real *dvv, void *user_data);
1421 status_t message_cb(message_t *msg, void *user_data);
1422 status_t interrupt_cb(integer *interrupt_status, void *user_data);
1423
1424 //=============================================================================
1425 /*!
1426  *
1427  */
1428 //=============================================================================
1429
1430 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
1431
1432   MESSAGE("BLSURFPlugin_BLSURF::Compute");
1433
1434   // Fix problem with locales
1435   Kernel_Utils::Localizer aLocalizer;
1436
1437   if ( !compute( aMesh, aShape ))
1438     return false;
1439
1440   if ( _haveViscousLayers )
1441   {
1442     // Compute viscous layers
1443
1444     TopTools_MapOfShape map;
1445     for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
1446     {
1447       const TopoDS_Face& f = TopoDS::Face(face_iter.Current());
1448       if ( !map.Add( f )) continue;
1449       SMESH_ProxyMesh::Ptr viscousMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, f );
1450       if ( !viscousMesh )
1451         return false; // error in StdMeshers_ViscousLayers2D::Compute()
1452
1453       // Compute BLSURF mesh on viscous layers
1454
1455       if ( viscousMesh->NbProxySubMeshes() > 0 )
1456         if ( !compute( aMesh, f, viscousMesh ))
1457           return false;
1458     }
1459
1460     // Re-compute BLSURF mesh on the rest faces if the mesh was cleared
1461
1462     for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
1463     {
1464       const TopoDS_Face& f = TopoDS::Face(face_iter.Current());
1465       SMESH_subMesh* fSM = aMesh.GetSubMesh( f );
1466       if ( fSM->IsMeshComputed() ) continue;
1467
1468       if ( !compute( aMesh, aShape ))
1469         return false;
1470       break;
1471     }
1472   }
1473   return true;
1474 }
1475
1476 //=============================================================================
1477 /*!
1478  *
1479  */
1480 //=============================================================================
1481
1482 bool BLSURFPlugin_BLSURF::compute(SMESH_Mesh&          aMesh,
1483                                   const TopoDS_Shape&  aShape,
1484                                   SMESH_ProxyMesh::Ptr viscousMesh)
1485 {
1486   /* create a distene context (generic object) */
1487   status_t status = STATUS_ERROR;
1488
1489   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1490   SMESH_MesherHelper helper( aMesh );
1491   // do not call helper.IsQuadraticSubMesh() because sub-meshes
1492   // may be cleaned and helper.myTLinkNodeMap gets invalid in such a case
1493   bool haveQuadraticSubMesh = SMESH_MesherHelper( aMesh ).IsQuadraticSubMesh( aShape );
1494   bool quadraticSubMeshAndViscousLayer = false;
1495   bool needMerge = false;
1496   typedef set< SMESHDS_SubMesh*, ShapeTypeCompare > TSubMeshSet;
1497   TSubMeshSet edgeSubmeshes, proxySubmeshes;
1498   TSubMeshSet& mergeSubmeshes = edgeSubmeshes;
1499
1500   TopTools_IndexedMapOfShape fmap;
1501   TopTools_IndexedMapOfShape emap;
1502   TopTools_IndexedMapOfShape pmap;
1503
1504   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1505 #ifndef WNT
1506   feclearexcept( FE_ALL_EXCEPT );
1507   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
1508 #endif
1509
1510   context_t *ctx =  context_new();
1511
1512   /* Set the message callback in the working context */
1513   context_set_message_callback(ctx, message_cb, &_comment);
1514 #ifdef WITH_SMESH_CANCEL_COMPUTE
1515   _compute_canceled = false;
1516   context_set_interrupt_callback(ctx, interrupt_cb, this);
1517 #else
1518   context_set_interrupt_callback(ctx, interrupt_cb, NULL);
1519 #endif
1520
1521   /* create the CAD object we will work on. It is associated to the context ctx. */
1522   cad_t *c     = cad_new(ctx);
1523   dcad_t *dcad = dcad_new(c);
1524
1525   FacesWithSizeMap.Clear();
1526   FaceId2SizeMap.clear();
1527   FaceId2ClassAttractor.clear();
1528   FaceIndex2ClassAttractor.clear();
1529   EdgesWithSizeMap.Clear();
1530   EdgeId2SizeMap.clear();
1531   VerticesWithSizeMap.Clear();
1532   VertexId2SizeMap.clear();
1533
1534   /* Now fill the CAD object with data from your CAD
1535    * environement. This is the most complex part of a successfull
1536    * integration.
1537    */
1538
1539   // PreCAD
1540   // If user requests it, send the CAD through Distene preprocessor : PreCAD
1541   cad_t *cleanc = NULL; // preprocessed cad
1542   precad_session_t *pcs = precad_session_new(ctx);
1543   precad_data_set_cad(pcs, c);
1544
1545   cadsurf_session_t *css = cadsurf_session_new(ctx);
1546
1547   // an object that correctly deletes all cadsurf objects at destruction
1548   BLSURF_Cleaner cleaner( ctx,css,c,dcad );
1549
1550   MESSAGE("BEGIN SetParameters");
1551   bool use_precad = false;
1552   SetParameters(
1553                 // #if BLSURF_VERSION_LONG >= "3.1.1"
1554                 //     c,
1555                 // #endif
1556                 _hypothesis, css, pcs, aMesh, &use_precad);
1557   MESSAGE("END SetParameters");
1558
1559   haveQuadraticSubMesh = haveQuadraticSubMesh || (_hypothesis != NULL && _hypothesis->GetQuadraticMesh());
1560   helper.SetIsQuadratic( haveQuadraticSubMesh );
1561
1562   // To remove as soon as quadratic mesh is allowed - BEGIN
1563   // GDD: Viscous layer is not allowed with quadratic mesh
1564   if (_haveViscousLayers && haveQuadraticSubMesh ) {
1565     quadraticSubMeshAndViscousLayer = true;
1566     _haveViscousLayers = !haveQuadraticSubMesh;
1567     _comment += "Warning: Viscous layer is not possible with a quadratic mesh, it is ignored.";
1568     error(COMPERR_WARNING, _comment);
1569   }
1570   // To remove as soon as quadratic mesh is allowed - END
1571
1572   // needed to prevent the opencascade memory managmement from freeing things
1573   vector<Handle(Geom2d_Curve)> curves;
1574   vector<Handle(Geom_Surface)> surfaces;
1575
1576   fmap.Clear();
1577   emap.Clear();
1578   pmap.Clear();
1579   FaceId2PythonSmp.clear();
1580   EdgeId2PythonSmp.clear();
1581   VertexId2PythonSmp.clear();
1582
1583   /****************************************************************************************
1584                                           FACES
1585   *****************************************************************************************/
1586   int iface = 0;
1587   string bad_end = "return";
1588   int faceKey = -1;
1589   TopTools_IndexedMapOfShape _map;
1590   TopExp::MapShapes(aShape,TopAbs_VERTEX,_map);
1591   int ienf = _map.Extent();
1592
1593   assert(Py_IsInitialized());
1594   PyGILState_STATE gstate;
1595   gstate = PyGILState_Ensure();
1596
1597   string theSizeMapStr;
1598
1599   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next())
1600   {
1601     TopoDS_Face f = TopoDS::Face(face_iter.Current());
1602
1603     // make INTERNAL face oriented FORWARD (issue 0020993)
1604     if (f.Orientation() != TopAbs_FORWARD && f.Orientation() != TopAbs_REVERSED )
1605       f.Orientation(TopAbs_FORWARD);
1606
1607     SMESH_subMesh* fSM = aMesh.GetSubMesh( f );
1608     if ( !viscousMesh && fSM->IsMeshComputed() )
1609       continue; // mesh on f already computed
1610
1611     if (fmap.FindIndex(f) > 0)
1612       continue;
1613     iface = fmap.Add(f);
1614
1615     surfaces.push_back(BRep_Tool::Surface(f));
1616
1617     /* create an object representing the face for cadsurf */
1618     /* where face_id is an integer identifying the face.
1619      * surf_function is the function that defines the surface
1620      * (For this face, it will be called by cadsurf with your_face_object_ptr
1621      * as last parameter.
1622      */
1623     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
1624
1625     /* by default a face has no tag (color).
1626        The following call sets it to the same value as the face_id : */
1627     cad_face_set_tag(fce, iface);
1628
1629     /* Set face orientation (optional if you want a well oriented output mesh)*/
1630     if(f.Orientation() != TopAbs_FORWARD)
1631       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
1632     else
1633       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
1634
1635     if (HasSizeMapOnFace && !use_precad)
1636     {
1637       // -----------------
1638       // Classic size map
1639       // -----------------
1640       faceKey = FacesWithSizeMap.FindIndex(f);
1641
1642
1643       if (FaceId2SizeMap.find(faceKey)!=FaceId2SizeMap.end()) {
1644         MESSAGE("A size map is defined on face :"<<faceKey)
1645           theSizeMapStr = FaceId2SizeMap[faceKey];
1646         // check if function ends with "return"
1647         if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1648           continue;
1649         // Expr To Python function, verification is performed at validation in GUI
1650         PyObject * obj = NULL;
1651         obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1652         Py_DECREF(obj);
1653         PyObject * func = NULL;
1654         func = PyObject_GetAttrString(main_mod, "f");
1655         FaceId2PythonSmp[iface]=func;
1656         FaceId2SizeMap.erase(faceKey);
1657       }
1658
1659       // Specific size map = Attractor
1660       std::map<int,std::vector<double> >::iterator attractor_iter = FaceId2AttractorCoords.begin();
1661
1662       for (; attractor_iter != FaceId2AttractorCoords.end(); ++attractor_iter) {
1663         if (attractor_iter->first == faceKey) {
1664           MESSAGE("Face indice: " << iface);
1665           MESSAGE("Adding attractor");
1666
1667           double xyzCoords[3]  = {attractor_iter->second[2],
1668                                   attractor_iter->second[3],
1669                                   attractor_iter->second[4]};
1670
1671           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
1672           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
1673           BRepClass_FaceClassifier scl(f,P,1e-7);
1674           // OCC 6.3sp6 : scl.Perform() is bugged. The function was rewritten
1675           // BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
1676           // OCC 6.5.2: scl.Perform() is not bugged anymore
1677           scl.Perform(f, P, 1e-7);
1678           TopAbs_State result = scl.State();
1679           MESSAGE("Position of point on face: "<<result);
1680           if ( result == TopAbs_OUT )
1681             MESSAGE("Point is out of face: node is not created");
1682           if ( result == TopAbs_UNKNOWN )
1683             MESSAGE("Point position on face is unknown: node is not created");
1684           if ( result == TopAbs_ON )
1685             MESSAGE("Point is on border of face: node is not created");
1686           if ( result == TopAbs_IN )
1687           {
1688             // Point is inside face and not on border
1689             MESSAGE("Point is in face: node is created");
1690             double uvCoords[2] = {attractor_iter->second[0],attractor_iter->second[1]};
1691             ienf++;
1692             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
1693             cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
1694             cad_point_set_tag(point_p, ienf);
1695           }
1696           FaceId2AttractorCoords.erase(faceKey);
1697         }
1698       }
1699
1700       // -----------------
1701       // Class Attractors
1702       // -----------------
1703       std::map<int,BLSURFPlugin_Attractor* >::iterator clAttractor_iter = FaceId2ClassAttractor.find(faceKey);
1704       if (clAttractor_iter != FaceId2ClassAttractor.end()){
1705         MESSAGE("Face indice: " << iface);
1706         MESSAGE("Adding attractor");
1707         FaceIndex2ClassAttractor[iface]=clAttractor_iter->second;
1708         FaceId2ClassAttractor.erase(clAttractor_iter);
1709       }
1710     } // if (HasSizeMapOnFace && !use_precad)
1711
1712       // ------------------
1713       // Enforced Vertices
1714       // ------------------
1715     faceKey = FacesWithEnforcedVertices.FindIndex(f);
1716     std::map<int,BLSURFPlugin_Hypothesis::TEnfVertexCoordsList >::const_iterator evmIt = FaceId2EnforcedVertexCoords.find(faceKey);
1717     if (evmIt != FaceId2EnforcedVertexCoords.end()) {
1718       MESSAGE("Some enforced vertices are defined");
1719       BLSURFPlugin_Hypothesis::TEnfVertexCoordsList evl;
1720       MESSAGE("Face indice: " << iface);
1721       MESSAGE("Adding enforced vertices");
1722       evl = evmIt->second;
1723       MESSAGE("Number of vertices to add: "<< evl.size());
1724       BLSURFPlugin_Hypothesis::TEnfVertexCoordsList::const_iterator evlIt = evl.begin();
1725       for (; evlIt != evl.end(); ++evlIt) {
1726         BLSURFPlugin_Hypothesis::TEnfVertexCoords xyzCoords;
1727         xyzCoords.push_back(evlIt->at(2));
1728         xyzCoords.push_back(evlIt->at(3));
1729         xyzCoords.push_back(evlIt->at(4));
1730         MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
1731         gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
1732         BRepClass_FaceClassifier scl(f,P,1e-7);
1733         // OCC 6.3sp6 : scl.Perform() is bugged. The function was rewritten
1734         // BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
1735         // OCC 6.5.2: scl.Perform() is not bugged anymore
1736         scl.Perform(f, P, 1e-7);
1737         TopAbs_State result = scl.State();
1738         MESSAGE("Position of point on face: "<<result);
1739         if ( result == TopAbs_OUT ) {
1740           MESSAGE("Point is out of face: node is not created");
1741           if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
1742             EnfVertexCoords2ProjVertex.erase(xyzCoords);
1743             EnfVertexCoords2EnfVertexList.erase(xyzCoords);
1744           }
1745         }
1746         if ( result == TopAbs_UNKNOWN ) {
1747           MESSAGE("Point position on face is unknown: node is not created");
1748           if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
1749             EnfVertexCoords2ProjVertex.erase(xyzCoords);
1750             EnfVertexCoords2EnfVertexList.erase(xyzCoords);
1751           }
1752         }
1753         if ( result == TopAbs_ON ) {
1754           MESSAGE("Point is on border of face: node is not created");
1755           if (EnfVertexCoords2ProjVertex.find(xyzCoords) != EnfVertexCoords2ProjVertex.end()) {
1756             EnfVertexCoords2ProjVertex.erase(xyzCoords);
1757             EnfVertexCoords2EnfVertexList.erase(xyzCoords);
1758           }
1759         }
1760         if ( result == TopAbs_IN )
1761         {
1762           // Point is inside face and not on border
1763           MESSAGE("Point is in face: node is created");
1764           double uvCoords[2]   = {evlIt->at(0),evlIt->at(1)};
1765           ienf++;
1766           MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
1767           cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
1768           int tag = 0;
1769           std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(xyzCoords);
1770           if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end() &&
1771               !enfCoordsIt->second.empty() )
1772           {
1773             // to merge nodes of an INTERNAL vertex belonging to several faces
1774             TopoDS_Vertex     v = (*enfCoordsIt->second.begin())->vertex;
1775             if ( v.IsNull() ) v = (*enfCoordsIt->second.rbegin())->vertex;
1776             if ( !v.IsNull() ) {
1777               tag = pmap.Add( v );
1778               SMESH_subMesh* vSM = aMesh.GetSubMesh( v );
1779               vSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1780               mergeSubmeshes.insert( vSM->GetSubMeshDS() );
1781               //if ( tag != pmap.Extent() )
1782               needMerge = true;
1783             }
1784           }
1785           if ( tag == 0 ) tag = ienf;
1786           cad_point_set_tag(point_p, tag);
1787         }
1788       }
1789       FaceId2EnforcedVertexCoords.erase(faceKey);
1790
1791     }
1792
1793     /****************************************************************************************
1794                                            EDGES
1795                         now create the edges associated to this face
1796     *****************************************************************************************/
1797     int edgeKey = -1;
1798     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next())
1799     {
1800       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
1801       int ic = emap.FindIndex(e);
1802       if (ic <= 0)
1803         ic = emap.Add(e);
1804
1805       const SMESH_ProxyMesh::SubMesh* proxySubMesh = 0;
1806       if ( viscousMesh )
1807         proxySubMesh = viscousMesh->GetProxySubMesh( e );
1808
1809       double tmin,tmax;
1810       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
1811       if ( proxySubMesh )
1812         curves.back() = new TProxyPCurve( proxySubMesh, curves.back() );
1813
1814       if (HasSizeMapOnEdge){
1815         edgeKey = EdgesWithSizeMap.FindIndex(e);
1816         if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end()) {
1817           MESSAGE("Sizemap defined on edge with index " << edgeKey);
1818           theSizeMapStr = EdgeId2SizeMap[edgeKey];
1819           if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1820             continue;
1821           // Expr To Python function, verification is performed at validation in GUI
1822           PyObject * obj = NULL;
1823           obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1824           Py_DECREF(obj);
1825           PyObject * func = NULL;
1826           func = PyObject_GetAttrString(main_mod, "f");
1827           EdgeId2PythonSmp[ic]=func;
1828           EdgeId2SizeMap.erase(edgeKey);
1829         }
1830       }
1831       /* data of nodes existing on the edge */
1832       StdMeshers_FaceSidePtr nodeData;
1833       SMESH_subMesh* sm = aMesh.GetSubMesh( e );
1834       if ( !sm->IsEmpty() )
1835       {
1836         SMESH_subMeshIteratorPtr subsmIt = sm->getDependsOnIterator( /*includeSelf=*/true,
1837                                                                      /*complexFirst=*/false);
1838         while ( subsmIt->more() )
1839           edgeSubmeshes.insert( subsmIt->next()->GetSubMeshDS() );
1840
1841         if ( proxySubMesh )
1842         { // protect proxySubmeshes from clearing
1843           proxySubmeshes.insert( (SMESHDS_SubMesh*) proxySubMesh );
1844           proxySubmeshes.insert( meshDS->MeshElements( e ));
1845           proxySubmeshes.insert( meshDS->MeshElements( helper.IthVertex( 0, e )));
1846           proxySubmeshes.insert( meshDS->MeshElements( helper.IthVertex( 1, e )));
1847         }
1848         nodeData.reset( new StdMeshers_FaceSide( f, e, &aMesh, /*isForwrd = */true,
1849                                                  /*ignoreMedium=*/haveQuadraticSubMesh,
1850                                                  viscousMesh));
1851         if ( nodeData->MissVertexNode() )
1852           return error(COMPERR_BAD_INPUT_MESH,"No node on vertex");
1853
1854         const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
1855         if ( !nodeDataVec.empty() )
1856         {
1857           if ( Abs( nodeDataVec[0].param - tmin ) > Abs( nodeDataVec.back().param - tmin ))
1858           {
1859             nodeData->Reverse();
1860             nodeData->GetUVPtStruct(); // nodeData recomputes nodeDataVec
1861           }
1862           // tmin and tmax can change in case of viscous layer on an adjacent edge
1863           tmin = nodeDataVec.front().param;
1864           tmax = nodeDataVec.back().param;
1865         }
1866         else
1867         {
1868           cout << "---------------- Invalid nodeData" << endl;
1869           nodeData.reset();
1870         }
1871       }
1872
1873       /* attach the edge to the current cadsurf face */
1874       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
1875
1876       /* by default an edge has no tag (color).
1877          The following call sets it to the same value as the edge_id : */
1878       cad_edge_set_tag(edg, ic);
1879
1880       /* by default, an edge does not necessalry appear in the resulting mesh,
1881          unless the following property is set :
1882       */
1883       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
1884
1885       /* by default an edge is a boundary edge */
1886       if (e.Orientation() == TopAbs_INTERNAL)
1887         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
1888
1889       // pass existing nodes of sub-meshes to BLSURF
1890       if ( nodeData )
1891       {
1892         const std::vector<UVPtStruct>& nodeDataVec = nodeData->GetUVPtStruct();
1893         const int                      nbNodes     = nodeDataVec.size();
1894
1895         dcad_edge_discretization_t *dedge;
1896         dcad_get_edge_discretization(dcad, edg, &dedge);
1897         dcad_edge_discretization_set_vertex_count( dedge, nbNodes );
1898
1899         //cout << endl << " EDGE " << ic << endl;
1900         //cout << "tmin = "<<tmin << ", tmax = "<< tmax << endl;
1901         for ( int iN = 0; iN < nbNodes; ++iN )
1902         {
1903           const UVPtStruct& nData = nodeDataVec[ iN ];
1904           double t                = nData.param;
1905           real uv[2]              = { nData.u, nData.v };
1906           SMESH_TNodeXYZ nXYZ( nData.node );
1907           // cout << "\tt = " << t << " uv = ( " << uv[0] << ","<< uv[1] << " ) ID " << nData.node->GetID()
1908           //      << " Curve UV: " << curves.back()->Value( t ).X()
1909           //      << ", "          << curves.back()->Value( t ).Y() << endl;
1910           dcad_edge_discretization_set_vertex_coordinates( dedge, iN+1, t, uv, nXYZ._xyz );
1911         }
1912         dcad_edge_discretization_set_property(dedge, DISTENE_DCAD_PROPERTY_REQUIRED);
1913       }
1914
1915       /****************************************************************************************
1916                                       VERTICES
1917       *****************************************************************************************/
1918
1919       int npts = 0;
1920       int ip1, ip2, *ip;
1921       gp_Pnt2d e0 = curves.back()->Value(tmin);
1922       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
1923       Standard_Real d1=0,d2=0;
1924
1925       int vertexKey = -1;
1926       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
1927         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
1928         ++npts;
1929         if (npts == 1){
1930           ip = &ip1;
1931           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1932         } else {
1933           ip = &ip2;
1934           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1935         }
1936         *ip = pmap.FindIndex(v);
1937         if(*ip <= 0)
1938           *ip = pmap.Add(v);
1939
1940         if (HasSizeMapOnVertex){
1941           vertexKey = VerticesWithSizeMap.FindIndex(v);
1942           if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
1943             theSizeMapStr = VertexId2SizeMap[vertexKey];
1944             //MESSAGE("VertexId2SizeMap[faceKey]: " << VertexId2SizeMap[vertexKey]);
1945             if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1946               continue;
1947             // Expr To Python function, verification is performed at validation in GUI
1948             PyObject * obj = NULL;
1949             obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1950             Py_DECREF(obj);
1951             PyObject * func = NULL;
1952             func = PyObject_GetAttrString(main_mod, "f");
1953             VertexId2PythonSmp[*ip]=func;
1954             VertexId2SizeMap.erase(vertexKey);   // do not erase if using a vector
1955           }
1956         }
1957       }
1958       if (npts != 2) {
1959         // should not happen
1960         MESSAGE("An edge does not have 2 extremities.");
1961       } else {
1962         if (d1 < d2) {
1963           // This defines the curves extremity connectivity
1964           cad_edge_set_extremities(edg, ip1, ip2);
1965           /* set the tag (color) to the same value as the extremity id : */
1966           cad_edge_set_extremities_tag(edg, ip1, ip2);
1967         }
1968         else {
1969           cad_edge_set_extremities(edg, ip2, ip1);
1970           cad_edge_set_extremities_tag(edg, ip2, ip1);
1971         }
1972       }
1973     } // for edge
1974   } //for face
1975
1976   // Clear mesh from already meshed edges if possible else
1977   // remember that merge is needed
1978   TSubMeshSet::iterator smIt = edgeSubmeshes.begin();
1979   for ( ; smIt != edgeSubmeshes.end(); ++smIt ) // loop on already meshed EDGEs
1980   {
1981     SMESHDS_SubMesh* smDS = *smIt;
1982     if ( !smDS ) continue;
1983     SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
1984     if ( nIt->more() )
1985     {
1986       const SMDS_MeshNode* n = nIt->next();
1987       if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
1988       {
1989         needMerge = true;
1990         // add existing medium nodes to helper
1991         if ( aMesh.NbEdges( ORDER_QUADRATIC ) > 0 )
1992         {
1993           SMDS_ElemIteratorPtr edgeIt = smDS->GetElements();
1994           while ( edgeIt->more() )
1995             helper.AddTLinks( static_cast<const SMDS_MeshEdge*>(edgeIt->next()));
1996         }
1997         continue;
1998       }
1999     }
2000     if ( proxySubmeshes.count( smDS ))
2001     {
2002       needMerge = true;
2003     }
2004     else
2005     {
2006       SMDS_ElemIteratorPtr eIt = smDS->GetElements();
2007       while ( eIt->more() ) meshDS->RemoveFreeElement( eIt->next(), smDS );
2008       SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2009       while ( nIt->more() ) meshDS->RemoveFreeNode( nIt->next(), smDS );
2010     }
2011   }
2012
2013
2014   if (use_precad) {
2015     /* Now launch the PreCAD process */
2016     status = precad_process(pcs);
2017     if(status != STATUS_OK){
2018       cout << "PreCAD processing failed with error code " << status << "\n";
2019     }
2020     else {
2021       // retrieve the pre-processed CAD object
2022       cleanc = precad_new_cad(pcs);
2023       if(!cleanc){
2024         cout << "Unable to retrieve PreCAD result \n";
2025       }
2026       cout << "PreCAD processing successfull \n";
2027
2028       // #if BLSURF_VERSION_LONG >= "3.1.1"
2029       //       /* We can now get the updated sizemaps (if any) */
2030       // //       if(geo_sizemap_e)
2031       // //         clean_geo_sizemap_e = precad_new_sizemap(pcs, geo_sizemap_e);
2032       // //
2033       // //       if(geo_sizemap_f)
2034       // //         clean_geo_sizemap_f = precad_new_sizemap(pcs, geo_sizemap_f);
2035       //
2036       //       if(iso_sizemap_p)
2037       //         clean_iso_sizemap_p = precad_new_sizemap(pcs, iso_sizemap_p);
2038       //
2039       //       if(iso_sizemap_e)
2040       //         clean_iso_sizemap_e = precad_new_sizemap(pcs, iso_sizemap_e);
2041       //
2042       //       if(iso_sizemap_f)
2043       //         clean_iso_sizemap_f = precad_new_sizemap(pcs, iso_sizemap_f);
2044       // #endif
2045     }
2046     // Now we can delete the PreCAD session
2047     precad_session_delete(pcs);
2048   }
2049
2050   cadsurf_data_set_dcad(css, dcad);
2051   if (cleanc) {
2052     // Give the pre-processed CAD object to the current BLSurf session
2053     cadsurf_data_set_cad(css, cleanc);
2054   }
2055   else {
2056     // Use the original one
2057     cadsurf_data_set_cad(css, c);
2058   }
2059
2060   std::cout << std::endl;
2061   std::cout << "Beginning of Surface Mesh generation" << std::endl;
2062   std::cout << std::endl;
2063
2064   try {
2065     OCC_CATCH_SIGNALS;
2066
2067     status = cadsurf_compute_mesh(css);
2068
2069   }
2070   catch ( std::exception& exc ) {
2071     _comment += exc.what();
2072   }
2073   catch (Standard_Failure& ex) {
2074     _comment += ex.DynamicType()->Name();
2075     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
2076       _comment += ": ";
2077       _comment += ex.GetMessageString();
2078     }
2079   }
2080   catch (...) {
2081     if ( _comment.empty() )
2082       _comment = "Exception in cadsurf_compute_mesh()";
2083   }
2084   if ( status != STATUS_OK) {
2085     // There was an error while meshing
2086     error(_comment);
2087   }
2088
2089   PyGILState_Release(gstate);
2090
2091   std::cout << std::endl;
2092   std::cout << "End of Surface Mesh generation" << std::endl;
2093   std::cout << std::endl;
2094
2095   mesh_t *msh = NULL;
2096   cadsurf_data_get_mesh(css, &msh);
2097   if(!msh){
2098     /* release the mesh object */
2099     cadsurf_data_regain_mesh(css, msh);
2100     return error(_comment);
2101   }
2102
2103   std::string GMFFileName = BLSURFPlugin_Hypothesis::GetDefaultGMFFile();
2104   if (_hypothesis)
2105     GMFFileName = _hypothesis->GetGMFFile();
2106   if (GMFFileName != "") {
2107     //     bool GMFFileMode = _hypothesis->GetGMFFileMode();
2108     bool asciiFound = (GMFFileName.find(".mesh",GMFFileName.length()-5) != std::string::npos);
2109     bool binaryFound = (GMFFileName.find(".meshb",GMFFileName.length()-6) != std::string::npos);
2110     if (!asciiFound && !binaryFound)
2111       GMFFileName.append(".mesh");
2112     mesh_write_mesh(msh, GMFFileName.c_str());
2113   }
2114
2115   /* retrieve mesh data (see meshgems/mesh.h) */
2116   integer nv, ne, nt, nq, vtx[4], tag;
2117   integer *evedg, *evtri, *evquad, type;
2118   real xyz[3];
2119
2120   mesh_get_vertex_count(msh, &nv);
2121   mesh_get_edge_count(msh, &ne);
2122   mesh_get_triangle_count(msh, &nt);
2123   mesh_get_quadrangle_count(msh, &nq);
2124
2125   evedg  = (integer *)mesh_calloc_generic_buffer(msh);
2126   evtri  = (integer *)mesh_calloc_generic_buffer(msh);
2127   evquad = (integer *)mesh_calloc_generic_buffer(msh);
2128
2129   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
2130   bool* tags = new bool[nv+1];
2131
2132   /* enumerated vertices */
2133   for(int iv=1;iv<=nv;iv++) {
2134     mesh_get_vertex_coordinates(msh, iv, xyz);
2135     mesh_get_vertex_tag(msh, iv, &tag);
2136     // Issue 0020656. Use vertex coordinates
2137     if ( tag > 0 && tag <= pmap.Extent() ) {
2138       TopoDS_Vertex v = TopoDS::Vertex(pmap(tag));
2139       double tol = BRep_Tool::Tolerance( v );
2140       gp_Pnt p = BRep_Tool::Pnt( v );
2141       if ( p.IsEqual( gp_Pnt( xyz[0], xyz[1], xyz[2]), 2*tol))
2142         xyz[0] = p.X(), xyz[1] = p.Y(), xyz[2] = p.Z();
2143       else
2144         tag = 0; // enforced or attracted vertex
2145     }
2146     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
2147
2148     // Create group of enforced vertices if requested
2149     BLSURFPlugin_Hypothesis::TEnfVertexCoords projVertex;
2150     projVertex.clear();
2151     projVertex.push_back((double)xyz[0]);
2152     projVertex.push_back((double)xyz[1]);
2153     projVertex.push_back((double)xyz[2]);
2154     std::map< BLSURFPlugin_Hypothesis::TEnfVertexCoords, BLSURFPlugin_Hypothesis::TEnfVertexList >::const_iterator enfCoordsIt = EnfVertexCoords2EnfVertexList.find(projVertex);
2155     if (enfCoordsIt != EnfVertexCoords2EnfVertexList.end()) {
2156       MESSAGE("Found enforced vertex @ " << xyz[0] << ", " << xyz[1] << ", " << xyz[2]);
2157       BLSURFPlugin_Hypothesis::TEnfVertexList::const_iterator enfListIt = enfCoordsIt->second.begin();
2158       BLSURFPlugin_Hypothesis::TEnfVertex *currentEnfVertex;
2159       for (; enfListIt != enfCoordsIt->second.end(); ++enfListIt) {
2160         currentEnfVertex = (*enfListIt);
2161         if (currentEnfVertex->grpName != "") {
2162           bool groupDone = false;
2163           SMESH_Mesh::GroupIteratorPtr grIt = aMesh.GetGroups();
2164           MESSAGE("currentEnfVertex->grpName: " << currentEnfVertex->grpName);
2165           MESSAGE("Parsing the groups of the mesh");
2166           while (grIt->more()) {
2167             SMESH_Group * group = grIt->next();
2168             if ( !group ) continue;
2169             MESSAGE("Group: " << group->GetName());
2170             SMESHDS_GroupBase* groupDS = group->GetGroupDS();
2171             if ( !groupDS ) continue;
2172             MESSAGE("group->SMDSGroup().GetType(): " << (groupDS->GetType()));
2173             MESSAGE("group->SMDSGroup().GetType()==SMDSAbs_Node: " << (groupDS->GetType()==SMDSAbs_Node));
2174             MESSAGE("currentEnfVertex->grpName.compare(group->GetStoreName())==0: " << (currentEnfVertex->grpName.compare(group->GetName())==0));
2175             if ( groupDS->GetType()==SMDSAbs_Node && currentEnfVertex->grpName.compare(group->GetName())==0) {
2176               SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
2177               aGroupDS->SMDSGroup().Add(nodes[iv]);
2178               MESSAGE("Node ID: " << nodes[iv]->GetID());
2179               // How can I inform the hypothesis ?
2180               //                 _hypothesis->AddEnfVertexNodeID(currentEnfVertex->grpName,nodes[iv]->GetID());
2181               groupDone = true;
2182               MESSAGE("Successfully added enforced vertex to existing group " << currentEnfVertex->grpName);
2183               break;
2184             }
2185           }
2186           if (!groupDone)
2187           {
2188             int groupId;
2189             SMESH_Group* aGroup = aMesh.AddGroup(SMDSAbs_Node, currentEnfVertex->grpName.c_str(), groupId);
2190             aGroup->SetName( currentEnfVertex->grpName.c_str() );
2191             SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
2192             aGroupDS->SMDSGroup().Add(nodes[iv]);
2193             MESSAGE("Successfully created enforced vertex group " << currentEnfVertex->grpName);
2194             groupDone = true;
2195           }
2196           if (!groupDone)
2197             throw SALOME_Exception(LOCALIZED("An enforced vertex node was not added to a group"));
2198         }
2199         else
2200           MESSAGE("Group name is empty: '"<<currentEnfVertex->grpName<<"' => group is not created");
2201       }
2202     }
2203
2204     // internal points are tagged to zero
2205     if(tag > 0 && tag <= pmap.Extent() ){
2206       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
2207       tags[iv] = false;
2208     } else {
2209       tags[iv] = true;
2210     }
2211   }
2212
2213   /* enumerate edges */
2214   for(int it=1;it<=ne;it++) {
2215     SMDS_MeshEdge* edg;
2216     mesh_get_edge_vertices(msh, it, vtx);
2217     mesh_get_edge_extra_vertices(msh, it, &type, evedg);
2218     mesh_get_edge_tag(msh, it, &tag);
2219     if (tags[vtx[0]]) {
2220       Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
2221       tags[vtx[0]] = false;
2222     };
2223     if (tags[vtx[1]]) {
2224       Set_NodeOnEdge(meshDS, nodes[vtx[1]], emap(tag));
2225       tags[vtx[1]] = false;
2226     };
2227     if (type == MESHGEMS_MESH_ELEMENT_TYPE_EDGE3) {
2228       // QUADRATIC EDGE
2229       if (tags[evedg[0]]) {
2230         Set_NodeOnEdge(meshDS, nodes[evedg[0]], emap(tag));
2231         tags[evedg[0]] = false;
2232       }
2233       edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]], nodes[evedg[0]]);
2234     }
2235     else {
2236       edg = helper.AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
2237     }
2238     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
2239   }
2240
2241   /* enumerate triangles */
2242   for(int it=1;it<=nt;it++) {
2243     SMDS_MeshFace* tri;
2244     mesh_get_triangle_vertices(msh, it, vtx);
2245     mesh_get_triangle_extra_vertices(msh, it, &type, evtri);
2246     mesh_get_triangle_tag(msh, it, &tag);
2247     if (tags[vtx[0]]) {
2248       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
2249       tags[vtx[0]] = false;
2250     };
2251     if (tags[vtx[1]]) {
2252       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
2253       tags[vtx[1]] = false;
2254     };
2255     if (tags[vtx[2]]) {
2256       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
2257       tags[vtx[2]] = false;
2258     };
2259     if (type == MESHGEMS_MESH_ELEMENT_TYPE_TRIA6) {
2260       // QUADRATIC TRIANGLE
2261       if (tags[evtri[0]]) {
2262         meshDS->SetNodeOnFace(nodes[evtri[0]], TopoDS::Face(fmap(tag)));
2263         tags[evtri[0]] = false;
2264       }
2265       if (tags[evtri[1]]) {
2266         meshDS->SetNodeOnFace(nodes[evtri[1]], TopoDS::Face(fmap(tag)));
2267         tags[evtri[1]] = false;
2268       }
2269       if (tags[evtri[2]]) {
2270         meshDS->SetNodeOnFace(nodes[evtri[2]], TopoDS::Face(fmap(tag)));
2271         tags[evtri[2]] = false;
2272       }
2273       tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]],
2274                             nodes[evtri[0]], nodes[evtri[1]], nodes[evtri[2]]);
2275     }
2276     else {
2277       tri = helper.AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
2278     }
2279     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
2280   }
2281
2282   /* enumerate quadrangles */
2283   for(int it=1;it<=nq;it++) {
2284     SMDS_MeshFace* quad;
2285     mesh_get_quadrangle_vertices(msh, it, vtx);
2286     mesh_get_quadrangle_extra_vertices(msh, it, &type, evquad);
2287     mesh_get_quadrangle_tag(msh, it, &tag);
2288     if (tags[vtx[0]]) {
2289       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
2290       tags[vtx[0]] = false;
2291     };
2292     if (tags[vtx[1]]) {
2293       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
2294       tags[vtx[1]] = false;
2295     };
2296     if (tags[vtx[2]]) {
2297       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
2298       tags[vtx[2]] = false;
2299     };
2300     if (tags[vtx[3]]) {
2301       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
2302       tags[vtx[3]] = false;
2303     };
2304     if (type == MESHGEMS_MESH_ELEMENT_TYPE_QUAD9) {
2305       // QUADRATIC QUADRANGLE
2306       std::cout << "This is a quadratic quadrangle" << std::endl;
2307       if (tags[evquad[0]]) {
2308         meshDS->SetNodeOnFace(nodes[evquad[0]], TopoDS::Face(fmap(tag)));
2309         tags[evquad[0]] = false;
2310       }
2311       if (tags[evquad[1]]) {
2312         meshDS->SetNodeOnFace(nodes[evquad[1]], TopoDS::Face(fmap(tag)));
2313         tags[evquad[1]] = false;
2314       }
2315       if (tags[evquad[2]]) {
2316         meshDS->SetNodeOnFace(nodes[evquad[2]], TopoDS::Face(fmap(tag)));
2317         tags[evquad[2]] = false;
2318       }
2319       if (tags[evquad[3]]) {
2320         meshDS->SetNodeOnFace(nodes[evquad[3]], TopoDS::Face(fmap(tag)));
2321         tags[evquad[3]] = false;
2322       }
2323       if (tags[evquad[4]]) {
2324         meshDS->SetNodeOnFace(nodes[evquad[4]], TopoDS::Face(fmap(tag)));
2325         tags[evquad[4]] = false;
2326       }
2327       quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]],
2328                              nodes[evquad[0]], nodes[evquad[1]], nodes[evquad[2]], nodes[evquad[3]],
2329                              nodes[evquad[4]]);
2330     }
2331     else {
2332       quad = helper.AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
2333     }
2334     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
2335   }
2336
2337   /* release the mesh object, the rest is released by cleaner */
2338   cadsurf_data_regain_mesh(css, msh);
2339
2340   delete [] nodes;
2341   delete [] tags;
2342
2343   if ( needMerge ) // sew mesh computed by BLSURF with pre-existing mesh
2344   {
2345     SMESH_MeshEditor editor( &aMesh );
2346     SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
2347     TIDSortedElemSet segementsOnEdge;
2348     TIDSortedNodeSet nodesOnEdge;
2349     TSubMeshSet::iterator smIt;
2350     SMESHDS_SubMesh* smDS;
2351     typedef SMDS_StdIterator< const SMDS_MeshNode*, SMDS_NodeIteratorPtr > TNodeIterator;
2352     double tol;
2353
2354     if ( !proxySubmeshes.empty() )
2355     {
2356       // merge proxy nodes with ones computed by BLSURF
2357       for ( smIt = proxySubmeshes.begin(); smIt != proxySubmeshes.end(); ++smIt )
2358       {
2359         if (! (smDS = *smIt) ) continue;
2360         nodesOnEdge.clear();
2361         nodesOnEdge.insert( TNodeIterator( smDS->GetNodes() ), TNodeIterator() );
2362         if (( smDS = meshDS->MeshElements( smDS->GetID() )))
2363         {
2364           nodesOnEdge.insert( TNodeIterator( smDS->GetNodes() ), TNodeIterator() );
2365           if ( nodesOnEdge.size() > 1 )
2366           {
2367             tol = SMESH_TNodeXYZ( *nodesOnEdge.begin() ).Distance( *nodesOnEdge.rbegin() );
2368             tol /= 50. * nodesOnEdge.size();
2369             editor.FindCoincidentNodes( nodesOnEdge, tol, nodeGroupsToMerge );
2370             editor.MergeNodes( nodeGroupsToMerge );
2371           }
2372         }
2373       }
2374     }
2375     // merge nodes on EDGE's with ones computed by BLSURF
2376     for ( smIt = mergeSubmeshes.begin(); smIt != mergeSubmeshes.end(); ++smIt )
2377     {
2378       if (! (smDS = *smIt) ) continue;
2379       getNodeGroupsToMerge( smDS, meshDS->IndexToShape((*smIt)->GetID()), nodeGroupsToMerge );
2380
2381       SMDS_ElemIteratorPtr segIt = smDS->GetElements();
2382       while ( segIt->more() )
2383         segementsOnEdge.insert( segIt->next() );
2384     }
2385     // merge nodes
2386     editor.MergeNodes( nodeGroupsToMerge );
2387
2388     // merge segments
2389     SMESH_MeshEditor::TListOfListOfElementsID equalSegments;
2390     editor.FindEqualElements( segementsOnEdge, equalSegments );
2391     editor.MergeElements( equalSegments );
2392
2393     // remove excess segments created on the boundary of viscous layers
2394     const SMDS_TypeOfPosition onFace = SMDS_TOP_FACE;
2395     for ( int i = 1; i <= emap.Extent(); ++i )
2396     {
2397       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( emap( i )))
2398       {
2399         SMDS_ElemIteratorPtr segIt = smDS->GetElements();
2400         while ( segIt->more() )
2401         {
2402           const SMDS_MeshElement* seg = segIt->next();
2403           if ( seg->GetNode(0)->GetPosition()->GetTypeOfPosition() == onFace ||
2404                seg->GetNode(1)->GetPosition()->GetTypeOfPosition() == onFace )
2405             meshDS->RemoveFreeElement( seg, smDS );
2406         }
2407       }
2408     }
2409   }
2410
2411   // SetIsAlwaysComputed( true ) to sub-meshes of EDGEs w/o mesh
2412   TopLoc_Location loc; double f,l;
2413   for (int i = 1; i <= emap.Extent(); i++)
2414     if ( SMESH_subMesh* sm = aMesh.GetSubMeshContaining( emap( i )))
2415       sm->SetIsAlwaysComputed( true );
2416   for (int i = 1; i <= pmap.Extent(); i++)
2417     if ( SMESH_subMesh* sm = aMesh.GetSubMeshContaining( pmap( i )))
2418       if ( !sm->IsMeshComputed() )
2419         sm->SetIsAlwaysComputed( true );
2420
2421   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
2422 #ifndef WNT
2423   if ( oldFEFlags > 0 )
2424     feenableexcept( oldFEFlags );
2425   feclearexcept( FE_ALL_EXCEPT );
2426 #endif
2427
2428   /*
2429     std::cout << "FacesWithSizeMap" << std::endl;
2430     FacesWithSizeMap.Statistics(std::cout);
2431     std::cout << "EdgesWithSizeMap" << std::endl;
2432     EdgesWithSizeMap.Statistics(std::cout);
2433     std::cout << "VerticesWithSizeMap" << std::endl;
2434     VerticesWithSizeMap.Statistics(std::cout);
2435     std::cout << "FacesWithEnforcedVertices" << std::endl;
2436     FacesWithEnforcedVertices.Statistics(std::cout);
2437   */
2438
2439   MESSAGE("END OF BLSURFPlugin_BLSURF::Compute()");
2440   return ( status == STATUS_OK && !quadraticSubMeshAndViscousLayer );
2441 }
2442
2443 //================================================================================
2444 /*!
2445  * \brief Terminates computation
2446  */
2447 //================================================================================
2448
2449 #ifdef WITH_SMESH_CANCEL_COMPUTE
2450 void BLSURFPlugin_BLSURF::CancelCompute()
2451 {
2452   _compute_canceled = true;
2453 }
2454 #endif
2455
2456 //=============================================================================
2457 /*!
2458  *  SetNodeOnEdge
2459  */
2460 //=============================================================================
2461
2462 void BLSURFPlugin_BLSURF::Set_NodeOnEdge(SMESHDS_Mesh* meshDS, SMDS_MeshNode* node, const TopoDS_Shape& ed) {
2463   const TopoDS_Edge edge = TopoDS::Edge(ed);
2464
2465   gp_Pnt pnt(node->X(), node->Y(), node->Z());
2466
2467   Standard_Real p0 = 0.0;
2468   Standard_Real p1 = 1.0;
2469   TopLoc_Location loc;
2470   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, loc, p0, p1);
2471
2472   if ( !loc.IsIdentity() ) pnt.Transform( loc.Transformation().Inverted() );
2473   GeomAPI_ProjectPointOnCurve proj(pnt, curve, p0, p1);
2474
2475   double pa = 0.;
2476   if ( proj.NbPoints() > 0 )
2477   {
2478     pa = (double)proj.LowerDistanceParameter();
2479     // Issue 0020656. Move node if it is too far from edge
2480     gp_Pnt curve_pnt = curve->Value( pa );
2481     double dist2     = pnt.SquareDistance( curve_pnt );
2482     double tol       = BRep_Tool::Tolerance( edge );
2483     if ( 1e-14 < dist2 && dist2 <= 1000*tol ) // large enough and within tolerance
2484     {
2485       curve_pnt.Transform( loc );
2486       meshDS->MoveNode( node, curve_pnt.X(), curve_pnt.Y(), curve_pnt.Z() );
2487     }
2488   }
2489 //   GProp_GProps LProps;
2490 //   BRepGProp::LinearProperties(ed, LProps);
2491 //   double lg = (double)LProps.Mass();
2492
2493   meshDS->SetNodeOnEdge(node, edge, pa);
2494 }
2495
2496 /* Curve definition function See cad_curv_t in file meshgems/cad.h for
2497  * more information.
2498  * NOTE : if when your CAD systems evaluates second
2499  * order derivatives it also computes first order derivatives and
2500  * function evaluation, you can optimize this example by making only
2501  * one CAD call and filling the necessary uv, dt, dtt arrays.
2502  */
2503 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
2504 {
2505   /* t is given. It contains the t (time) 1D parametric coordintaes
2506      of the point PreCAD/BLSurf is querying on the curve */
2507
2508   /* user_data identifies the edge PreCAD/BLSurf is querying
2509    * (see cad_edge_new later in this example) */
2510   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
2511
2512   if (uv){
2513    /* BLSurf is querying the function evaluation */
2514     gp_Pnt2d P;
2515     P=pargeo->Value(t);
2516     uv[0]=P.X(); uv[1]=P.Y();
2517   }
2518
2519   if(dt) {
2520    /* query for the first order derivatives */
2521     gp_Vec2d V1;
2522     V1=pargeo->DN(t,1);
2523     dt[0]=V1.X(); dt[1]=V1.Y();
2524   }
2525
2526   if(dtt){
2527     /* query for the second order derivatives */
2528     gp_Vec2d V2;
2529     V2=pargeo->DN(t,2);
2530     dtt[0]=V2.X(); dtt[1]=V2.Y();
2531   }
2532
2533   return STATUS_OK;
2534 }
2535
2536 /* Surface definition function.
2537  * See cad_surf_t in file meshgems/cad.h for more information.
2538  * NOTE : if when your CAD systems evaluates second order derivatives it also
2539  * computes first order derivatives and function evaluation, you can optimize
2540  * this example by making only one CAD call and filling the necessary xyz, du, dv, etc..
2541  * arrays.
2542  */
2543 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
2544                   real *duu, real *duv, real *dvv, void *user_data)
2545 {
2546   /* uv[2] is given. It contains the u,v coordinates of the point
2547    * PreCAD/BLSurf is querying on the surface */
2548
2549   /* user_data identifies the face PreCAD/BLSurf is querying (see
2550    * cad_face_new later in this example)*/
2551   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
2552
2553   if(xyz){
2554    gp_Pnt P;
2555    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
2556    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
2557   }
2558
2559   if(du && dv){
2560     gp_Pnt P;
2561     gp_Vec D1U,D1V;
2562
2563     geometry->D1(uv[0],uv[1],P,D1U,D1V);
2564     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
2565     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
2566   }
2567
2568   if(duu && duv && dvv){
2569
2570     gp_Pnt P;
2571     gp_Vec D1U,D1V;
2572     gp_Vec D2U,D2V,D2UV;
2573
2574     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
2575     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
2576     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
2577     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
2578   }
2579
2580   return STATUS_OK;
2581 }
2582
2583
2584 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
2585 {
2586   //MESSAGE("size_on_surface")
2587   if (FaceId2PythonSmp.count(face_id) != 0){
2588     //MESSAGE("A size map is used to calculate size on face : "<<face_id)
2589     PyObject * pyresult = NULL;
2590     PyObject* new_stderr = NULL;
2591     assert(Py_IsInitialized());
2592     PyGILState_STATE gstate;
2593     gstate = PyGILState_Ensure();
2594     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],(char*)"(f,f)",uv[0],uv[1]);
2595     real result;
2596     if ( pyresult != NULL) {
2597       result = PyFloat_AsDouble(pyresult);
2598       Py_DECREF(pyresult);
2599 //       *size = result;
2600     }
2601     else{
2602       fflush(stderr);
2603       string err_description="";
2604       new_stderr = newPyStdOut(err_description);
2605       PySys_SetObject((char*)"stderr", new_stderr);
2606       PyErr_Print();
2607       PySys_SetObject((char*)"stderr", PySys_GetObject((char*)"__stderr__"));
2608       Py_DECREF(new_stderr);
2609       MESSAGE("Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description);
2610       result = *((real*)user_data);
2611     }
2612     *size = result;
2613     PyGILState_Release(gstate);
2614   }
2615   else if (FaceIndex2ClassAttractor.count(face_id) !=0 && !FaceIndex2ClassAttractor[face_id]->Empty()){
2616 //    MESSAGE("attractor used on face :"<<face_id)
2617     // MESSAGE("List of attractor is not empty")
2618     // MESSAGE("Attractor empty : "<< FaceIndex2ClassAttractor[face_id]->Empty())
2619     real result = FaceIndex2ClassAttractor[face_id]->GetSize(uv[0],uv[1]);
2620     *size = result;
2621   }
2622   else {
2623     // MESSAGE("List of attractor is empty !!!")
2624     *size = *((real*)user_data);
2625   }
2626 //   std::cout << "Size_on_surface sur la face " << face_id << " donne une size de: " << *size << std::endl;
2627   return STATUS_OK;
2628 }
2629
2630 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
2631 {
2632   if (EdgeId2PythonSmp.count(edge_id) != 0){
2633     PyObject * pyresult = NULL;
2634     PyObject* new_stderr = NULL;
2635     assert(Py_IsInitialized());
2636     PyGILState_STATE gstate;
2637     gstate = PyGILState_Ensure();
2638     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],(char*)"(f)",t);
2639     real result;
2640     if ( pyresult != NULL) {
2641       result = PyFloat_AsDouble(pyresult);
2642       Py_DECREF(pyresult);
2643 //       *size = result;
2644     }
2645     else{
2646       fflush(stderr);
2647       string err_description="";
2648       new_stderr = newPyStdOut(err_description);
2649       PySys_SetObject((char*)"stderr", new_stderr);
2650       PyErr_Print();
2651       PySys_SetObject((char*)"stderr", PySys_GetObject((char*)"__stderr__"));
2652       Py_DECREF(new_stderr);
2653       MESSAGE("Can't evaluate f(" << t << ")" << " error is " << err_description);
2654       result = *((real*)user_data);
2655     }
2656     *size = result;
2657     PyGILState_Release(gstate);
2658   }
2659   else {
2660     *size = *((real*)user_data);
2661   }
2662   return STATUS_OK;
2663 }
2664
2665 status_t size_on_vertex(integer point_id, real *size, void *user_data)
2666 {
2667   if (VertexId2PythonSmp.count(point_id) != 0){
2668     PyObject * pyresult = NULL;
2669     PyObject* new_stderr = NULL;
2670     assert(Py_IsInitialized());
2671     PyGILState_STATE gstate;
2672     gstate = PyGILState_Ensure();
2673     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],(char*)"");
2674     real result;
2675     if ( pyresult != NULL) {
2676       result = PyFloat_AsDouble(pyresult);
2677       Py_DECREF(pyresult);
2678 //       *size = result;
2679     }
2680     else {
2681       fflush(stderr);
2682       string err_description="";
2683       new_stderr = newPyStdOut(err_description);
2684       PySys_SetObject((char*)"stderr", new_stderr);
2685       PyErr_Print();
2686       PySys_SetObject((char*)"stderr", PySys_GetObject((char*)"__stderr__"));
2687       Py_DECREF(new_stderr);
2688       MESSAGE("Can't evaluate f()" << " error is " << err_description);
2689       result = *((real*)user_data);
2690     }
2691     *size = result;
2692     PyGILState_Release(gstate);
2693   }
2694   else {
2695     *size = *((real*)user_data);
2696   }
2697  return STATUS_OK;
2698 }
2699
2700 /*
2701  * The following function will be called for PreCAD/BLSurf message
2702  * printing.  See context_set_message_callback (later in this
2703  * template) for how to set user_data.
2704  */
2705 status_t message_cb(message_t *msg, void *user_data)
2706 {
2707   integer errnumber = 0;
2708   char *desc;
2709   message_get_number(msg, &errnumber);
2710   message_get_description(msg, &desc);
2711   string err( desc );
2712   if ( errnumber < 0 || err.find("license") != string::npos ) {
2713     string * error = (string*)user_data;
2714 //   if ( !error->empty() )
2715 //     *error += "\n";
2716     // remove ^A from the tail
2717     int len = strlen( desc );
2718     while (len > 0 && desc[len-1] != '\n')
2719       len--;
2720     error->append( desc, len );
2721   }
2722   else {
2723       std::cout << desc << std::endl;
2724   }
2725   return STATUS_OK;
2726 }
2727
2728 /* This is the interrupt callback. PreCAD/BLSurf will call this
2729  * function regularily. See the file meshgems/interrupt.h
2730  */
2731 status_t interrupt_cb(integer *interrupt_status, void *user_data)
2732 {
2733   integer you_want_to_continue = 1;
2734 #ifdef WITH_SMESH_CANCEL_COMPUTE
2735   BLSURFPlugin_BLSURF* tmp = (BLSURFPlugin_BLSURF*)user_data;
2736   you_want_to_continue = !tmp->computeCanceled();
2737 #endif
2738
2739   if(you_want_to_continue)
2740   {
2741     *interrupt_status = INTERRUPT_CONTINUE;
2742     return STATUS_OK;
2743   }
2744   else /* you want to stop BLSurf */
2745   {
2746     *interrupt_status = INTERRUPT_STOP;
2747     return STATUS_ERROR;
2748   }
2749 }
2750
2751 //=============================================================================
2752 /*!
2753  *
2754  */
2755 //=============================================================================
2756 bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh&         aMesh,
2757                                    const TopoDS_Shape& aShape,
2758                                    MapShapeNbElems&    aResMap)
2759 {
2760   double diagonal       = aMesh.GetShapeDiagonalSize();
2761   double bbSegmentation = _gen->GetBoundaryBoxSegmentation();
2762   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
2763   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize(diagonal, bbSegmentation);
2764   bool   _phySizeRel    = BLSURFPlugin_Hypothesis::GetDefaultPhySizeRel();
2765   //int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
2766   double _angleMesh     = BLSURFPlugin_Hypothesis::GetDefaultAngleMesh();
2767   bool   _quadAllowed   = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
2768   if(_hypothesis) {
2769     _physicalMesh  = (int) _hypothesis->GetPhysicalMesh();
2770     _phySizeRel         = _hypothesis->IsPhySizeRel();
2771     if ( _hypothesis->GetPhySize() > 0)
2772       _phySize          = _phySizeRel ? diagonal*_hypothesis->GetPhySize() : _hypothesis->GetPhySize();
2773     //_geometricMesh = (int) hyp->GetGeometricMesh();
2774     if (_hypothesis->GetAngleMesh() > 0)
2775       _angleMesh        = _hypothesis->GetAngleMesh();
2776     _quadAllowed        = _hypothesis->GetQuadAllowed();
2777   } else {
2778     //0020968: EDF1545 SMESH: Problem in the creation of a mesh group on geometry
2779     // GetDefaultPhySize() sometimes leads to computation failure
2780     _phySize = aMesh.GetShapeDiagonalSize() / _gen->GetBoundaryBoxSegmentation();
2781     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
2782   }
2783
2784   bool IsQuadratic = _quadraticMesh;
2785
2786   // ----------------
2787   // evaluate 1D
2788   // ----------------
2789   TopTools_DataMapOfShapeInteger EdgesMap;
2790   double fullLen = 0.0;
2791   double fullNbSeg = 0;
2792   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
2793     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
2794     if( EdgesMap.IsBound(E) )
2795       continue;
2796     SMESH_subMesh *sm = aMesh.GetSubMesh(E);
2797     double aLen = SMESH_Algo::EdgeLength(E);
2798     fullLen += aLen;
2799     int nb1d = 0;
2800     if(_physicalMesh==1) {
2801        nb1d = (int)( aLen/_phySize + 1 );
2802     }
2803     else {
2804       // use geometry
2805       double f,l;
2806       Handle(Geom_Curve) C = BRep_Tool::Curve(E,f,l);
2807       double fullAng = 0.0;
2808       double dp = (l-f)/200;
2809       gp_Pnt P1,P2,P3;
2810       C->D0(f,P1);
2811       C->D0(f+dp,P2);
2812       gp_Vec V1(P1,P2);
2813       for(int j=2; j<=200; j++) {
2814         C->D0(f+dp*j,P3);
2815         gp_Vec V2(P2,P3);
2816         fullAng += fabs(V1.Angle(V2));
2817         V1 = V2;
2818         P2 = P3;
2819       }
2820       nb1d = (int)( fullAng/_angleMesh + 1 );
2821     }
2822     fullNbSeg += nb1d;
2823     std::vector<int> aVec(SMDSEntity_Last);
2824     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
2825     if( IsQuadratic > 0 ) {
2826       aVec[SMDSEntity_Node] = 2*nb1d - 1;
2827       aVec[SMDSEntity_Quad_Edge] = nb1d;
2828     }
2829     else {
2830       aVec[SMDSEntity_Node] = nb1d - 1;
2831       aVec[SMDSEntity_Edge] = nb1d;
2832     }
2833     aResMap.insert(std::make_pair(sm,aVec));
2834     EdgesMap.Bind(E,nb1d);
2835   }
2836   double ELen = fullLen/fullNbSeg;
2837   // ----------------
2838   // evaluate 2D
2839   // ----------------
2840   // try to evaluate as in MEFISTO
2841   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
2842     TopoDS_Face F = TopoDS::Face( exp.Current() );
2843     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
2844     GProp_GProps G;
2845     BRepGProp::SurfaceProperties(F,G);
2846     double anArea = G.Mass();
2847     int nb1d = 0;
2848     std::vector<int> nb1dVec;
2849     for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) {
2850       int nbSeg = EdgesMap.Find(exp1.Current());
2851       nb1d += nbSeg;
2852       nb1dVec.push_back( nbSeg );
2853     }
2854     int nbQuad = 0;
2855     int nbTria = (int) ( anArea/( ELen*ELen*sqrt(3.) / 4 ) );
2856     int nbNodes = (int) ( ( nbTria*3 - (nb1d-1)*2 ) / 6 + 1 );
2857     if ( _quadAllowed )
2858     {
2859       if ( nb1dVec.size() == 4 ) // quadrangle geom face
2860       {
2861         int n1 = nb1dVec[0], n2 = nb1dVec[ nb1dVec[1] == nb1dVec[0] ? 2 : 1 ];
2862         nbQuad = n1 * n2;
2863         nbNodes = (n1 + 1) * (n2 + 1);
2864         nbTria = 0;
2865       }
2866       else
2867       {
2868         nbTria = nbQuad = nbTria / 3 + 1;
2869       }
2870     }
2871     std::vector<int> aVec(SMDSEntity_Last,0);
2872     if( IsQuadratic ) {
2873       int nb1d_in = (nbTria*3 - nb1d) / 2;
2874       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
2875       aVec[SMDSEntity_Quad_Triangle] = nbTria;
2876       aVec[SMDSEntity_Quad_Quadrangle] = nbQuad;
2877     }
2878     else {
2879       aVec[SMDSEntity_Node] = nbNodes;
2880       aVec[SMDSEntity_Triangle] = nbTria;
2881       aVec[SMDSEntity_Quadrangle] = nbQuad;
2882     }
2883     aResMap.insert(std::make_pair(sm,aVec));
2884   }
2885
2886   // ----------------
2887   // evaluate 3D
2888   // ----------------
2889   GProp_GProps G;
2890   BRepGProp::VolumeProperties(aShape,G);
2891   double aVolume = G.Mass();
2892   double tetrVol = 0.1179*ELen*ELen*ELen;
2893   int nbVols  = int(aVolume/tetrVol);
2894   int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
2895   std::vector<int> aVec(SMDSEntity_Last);
2896   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
2897   if( IsQuadratic ) {
2898     aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
2899     aVec[SMDSEntity_Quad_Tetra] = nbVols;
2900   }
2901   else {
2902     aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
2903     aVec[SMDSEntity_Tetra] = nbVols;
2904   }
2905   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
2906   aResMap.insert(std::make_pair(sm,aVec));
2907
2908   return true;
2909 }