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