src/BLSURFPlugin/BLSURFPlugin_BLSURF.cxx

   1 //  Copyright (C) 2007-2008  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
  29 #include <structmember.h>
  30
  31
  32 #include <SMESH_Gen.hxx>
  33 #include <SMESH_Mesh.hxx>
  34 #include <SMESH_ControlsDef.hxx>
  35
  36 #include <SMESHDS_Mesh.hxx>
  37 #include <SMDS_MeshElement.hxx>
  38 #include <SMDS_MeshNode.hxx>
  39
  40 #include <utilities.h>
  41
  42 #include <limits>
  43 #include <list>
  44 #include <vector>
  45 #include <cstdlib>
  46
  47 #include <BRep_Tool.hxx>
  48 #include <TopExp.hxx>
  49 #include <TopExp_Explorer.hxx>
  50 #include <TopoDS.hxx>
  51 #include <NCollection_Map.hxx>
  52 #include <Standard_ErrorHandler.hxx>
  53
  54 extern "C"{
  55 #include <distene/api.h>
  56 }
  57
  58 #include <Geom_Surface.hxx>
  59 #include <Handle_Geom_Surface.hxx>
  60 #include <Geom2d_Curve.hxx>
  61 #include <Handle_Geom2d_Curve.hxx>
  62 #include <Geom_Curve.hxx>
  63 #include <Handle_Geom_Curve.hxx>
  64 #include <TopoDS_Vertex.hxx>
  65 #include <TopoDS_Edge.hxx>
  66 #include <TopoDS_Wire.hxx>
  67 #include <TopoDS_Face.hxx>
  68 #include <TopoDS_Shape.hxx>
  69 #include <gp_Pnt2d.hxx>
  70 #include <TopTools_IndexedMapOfShape.hxx>
  71 #include <BRepTools.hxx>
  72 #include <TopTools_DataMapOfShapeInteger.hxx>
  73 #include <GProp_GProps.hxx>
  74 #include <BRepGProp.hxx>
  75
  76 #ifndef WNT
  77 #include <fenv.h>
  78 #endif
  79
  80 #include <GeomAPI_ProjectPointOnCurve.hxx>
  81 #include <GeomAPI_ProjectPointOnSurf.hxx>
  82 #include <gp_XY.hxx>
  83 #include <gp_XYZ.hxx>
  84
  85 /* ==================================
  86  * ===========  PYTHON ==============
  87  * ==================================*/
  88
  89 typedef struct {
  90   PyObject_HEAD
  91   int softspace;
  92   std::string *out;
  93   } PyStdOut;
  94
  95 static void
  96 PyStdOut_dealloc(PyStdOut *self)
  97 {
  98   PyObject_Del(self);
  99 }
 100
 101 static PyObject *
 102 PyStdOut_write(PyStdOut *self, PyObject *args)
 103 {
 104   char *c;
 105   int l;
 106   if (!PyArg_ParseTuple(args, "t#:write",&c, &l))
 107     return NULL;
 108
 109   //std::cerr << c ;
 110   *(self->out)=*(self->out)+c;
 111
 112   Py_INCREF(Py_None);
 113   return Py_None;
 114 }
 115
 116 static PyMethodDef PyStdOut_methods[] = {
 117   {"write",  (PyCFunction)PyStdOut_write,  METH_VARARGS,
 118     PyDoc_STR("write(string) -> None")},
 119   {NULL,    NULL}   /* sentinel */
 120 };
 121
 122 static PyMemberDef PyStdOut_memberlist[] = {
 123   {"softspace", T_INT,  offsetof(PyStdOut, softspace), 0,
 124    "flag indicating that a space needs to be printed; used by print"},
 125   {NULL} /* Sentinel */
 126 };
 127
 128 static PyTypeObject PyStdOut_Type = {
 129   /* The ob_type field must be initialized in the module init function
 130    * to be portable to Windows without using C++. */
 131   PyObject_HEAD_INIT(NULL)
 132   0,      /*ob_size*/
 133   "PyOut",   /*tp_name*/
 134   sizeof(PyStdOut),  /*tp_basicsize*/
 135   0,      /*tp_itemsize*/
 136   /* methods */
 137   (destructor)PyStdOut_dealloc, /*tp_dealloc*/
 138   0,      /*tp_print*/
 139   0, /*tp_getattr*/
 140   0, /*tp_setattr*/
 141   0,      /*tp_compare*/
 142   0,      /*tp_repr*/
 143   0,      /*tp_as_number*/
 144   0,      /*tp_as_sequence*/
 145   0,      /*tp_as_mapping*/
 146   0,      /*tp_hash*/
 147         0,                      /*tp_call*/
 148         0,                      /*tp_str*/
 149         PyObject_GenericGetAttr,                      /*tp_getattro*/
 150         /* softspace is writable:  we must supply tp_setattro */
 151         PyObject_GenericSetAttr,    /* tp_setattro */
 152         0,                      /*tp_as_buffer*/
 153         Py_TPFLAGS_DEFAULT,     /*tp_flags*/
 154         0,                      /*tp_doc*/
 155         0,                      /*tp_traverse*/
 156         0,                      /*tp_clear*/
 157         0,                      /*tp_richcompare*/
 158         0,                      /*tp_weaklistoffset*/
 159         0,                      /*tp_iter*/
 160         0,                      /*tp_iternext*/
 161         PyStdOut_methods,                      /*tp_methods*/
 162         PyStdOut_memberlist,                      /*tp_members*/
 163         0,                      /*tp_getset*/
 164         0,                      /*tp_base*/
 165         0,                      /*tp_dict*/
 166         0,                      /*tp_descr_get*/
 167         0,                      /*tp_descr_set*/
 168         0,                      /*tp_dictoffset*/
 169         0,                      /*tp_init*/
 170         0,                      /*tp_alloc*/
 171         0,                      /*tp_new*/
 172         0,                      /*tp_free*/
 173         0,                      /*tp_is_gc*/
 174 };
 175
 176 PyObject * newPyStdOut( std::string& out )
 177 {
 178   PyStdOut *self;
 179   self = PyObject_New(PyStdOut, &PyStdOut_Type);
 180   if (self == NULL)
 181     return NULL;
 182   self->softspace = 0;
 183   self->out=&out;
 184   return (PyObject*)self;
 185 }
 186
 187
 188 ////////////////////////END PYTHON///////////////////////////
 189
 190 //////////////////MY MAPS////////////////////////////////////////
 191 std::map<int,string> FaceId2SizeMap;
 192 std::map<int,string> EdgeId2SizeMap;
 193 std::map<int,string> VertexId2SizeMap;
 194 std::map<int,PyObject*> FaceId2PythonSmp;
 195 std::map<int,PyObject*> EdgeId2PythonSmp;
 196 std::map<int,PyObject*> VertexId2PythonSmp;
 197
 198
 199 bool HasSizeMapOnFace=false;
 200 bool HasSizeMapOnEdge=false;
 201 bool HasSizeMapOnVertex=false;
 202
 203 //=============================================================================
 204 /*!
 205  *
 206  */
 207 //=============================================================================
 208
 209 BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF(int hypId, int studyId,
 210                                                SMESH_Gen* gen)
 211   : SMESH_2D_Algo(hypId, studyId, gen)
 212 {
 213   MESSAGE("BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF");
 214
 215   _name = "BLSURF";
 216   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
 217   _compatibleHypothesis.push_back("BLSURF_Parameters");
 218   _requireDescretBoundary = false;
 219   _onlyUnaryInput = false;
 220   _hypothesis = NULL;
 221
 222   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 223   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 224   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 225
 226   MESSAGE("studyid = " << _studyId);
 227
 228   myStudy = NULL;
 229   myStudy = aStudyMgr->GetStudyByID(_studyId);
 230   MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
 231
 232   /* Initialize the Python interpreter */
 233   assert(Py_IsInitialized());
 234   PyGILState_STATE gstate;
 235   gstate = PyGILState_Ensure();
 236
 237   main_mod = NULL;
 238   main_mod = PyImport_AddModule("__main__");
 239
 240   main_dict = NULL;
 241   main_dict = PyModule_GetDict(main_mod);
 242
 243   PyRun_SimpleString("from math import *");
 244   PyGILState_Release(gstate);
 245
 246   FaceId2SizeMap.clear();
 247   EdgeId2SizeMap.clear();
 248   VertexId2SizeMap.clear();
 249   FaceId2PythonSmp.clear();
 250   EdgeId2PythonSmp.clear();
 251   VertexId2PythonSmp.clear();
 252 }
 253
 254 //=============================================================================
 255 /*!
 256  *
 257  */
 258 //=============================================================================
 259
 260 BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF()
 261 {
 262   MESSAGE("BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF");
 263 }
 264
 265
 266 //=============================================================================
 267 /*!
 268  *
 269  */
 270 //=============================================================================
 271
 272 bool BLSURFPlugin_BLSURF::CheckHypothesis
 273                          (SMESH_Mesh&                          aMesh,
 274                           const TopoDS_Shape&                  aShape,
 275                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 276 {
 277   _hypothesis = NULL;
 278
 279   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 280   const SMESHDS_Hypothesis* theHyp;
 281
 282   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 283   int nbHyp = hyps.size();
 284   if (!nbHyp)
 285   {
 286     aStatus = SMESH_Hypothesis::HYP_OK;
 287     return true;  // can work with no hypothesis
 288   }
 289
 290   itl = hyps.begin();
 291   theHyp = (*itl); // use only the first hypothesis
 292
 293   string hypName = theHyp->GetName();
 294
 295   if (hypName == "BLSURF_Parameters")
 296   {
 297     _hypothesis = static_cast<const BLSURFPlugin_Hypothesis*> (theHyp);
 298     ASSERT(_hypothesis);
 299     if ( _hypothesis->GetPhysicalMesh() == BLSURFPlugin_Hypothesis::DefaultSize &&
 300          _hypothesis->GetGeometricMesh() == BLSURFPlugin_Hypothesis::DefaultGeom )
 301       //  hphy_flag = 0 and hgeo_flag = 0 is not allowed (spec)
 302       aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
 303     else
 304       aStatus = SMESH_Hypothesis::HYP_OK;
 305   }
 306   else
 307     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 308
 309   return aStatus == SMESH_Hypothesis::HYP_OK;
 310 }
 311
 312 //=============================================================================
 313 /*!
 314  * Pass parameters to BLSURF
 315  */
 316 //=============================================================================
 317
 318 inline std::string to_string(double d)
 319 {
 320    std::ostringstream o;
 321    o << d;
 322    return o.str();
 323 }
 324
 325 inline std::string to_string(int i)
 326 {
 327    std::ostringstream o;
 328    o << i;
 329    return o.str();
 330 }
 331
 332 double _smp_phy_size;
 333 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data);
 334 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data);
 335 status_t size_on_vertex(integer vertex_id, real *size, void *user_data);
 336
 337 double my_u_min=1e6,my_v_min=1e6,my_u_max=-1e6,my_v_max=-1e6;
 338
 339 /////////////////////////////////////////////////////////
 340 gp_XY getUV(const TopoDS_Face& face, const gp_XYZ& point)
 341 {
 342   Handle(Geom_Surface) surface = BRep_Tool::Surface(face);
 343   GeomAPI_ProjectPointOnSurf projector( point, surface );
 344   if ( !projector.IsDone() || projector.NbPoints()==0 )
 345     throw "Can't project";
 346
 347   Quantity_Parameter u,v;
 348   projector.LowerDistanceParameters(u,v);
 349   return gp_XY(u,v);
 350 }
 351 /////////////////////////////////////////////////////////
 352
 353 /////////////////////////////////////////////////////////
 354 double getT(const TopoDS_Edge& edge, const gp_XYZ& point)
 355 {
 356   Standard_Real f,l;
 357   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, f,l);
 358   GeomAPI_ProjectPointOnCurve projector( point, curve);
 359   if ( projector.NbPoints() == 0 )
 360     throw;
 361   return projector.LowerDistanceParameter();
 362 }
 363
 364 /////////////////////////////////////////////////////////
 365 TopoDS_Shape BLSURFPlugin_BLSURF::entryToShape(std::string entry)
 366 {
 367     MESSAGE("BLSURFPlugin_BLSURF::entryToShape"<<entry );
 368     TopoDS_Shape S = TopoDS_Shape();
 369     SALOMEDS::SObject_var aSO = myStudy->FindObjectID(entry.c_str());
 370     SALOMEDS::GenericAttribute_var anAttr;
 371     if (!aSO->_is_nil()){
 372       SALOMEDS::SObject_var aRefSObj;
 373       GEOM::GEOM_Object_var aShape;
 374       SALOMEDS::AttributeIOR_var myAttribute;
 375       CORBA::String_var myAttrValue;
 376       CORBA::Object_var myCorbaObj;
 377       // If selected object is a reference
 378       if ( aSO->ReferencedObject( aRefSObj ))
 379         aSO = aRefSObj;
 380       SALOMEDS::SComponent_var myFatherCpnt = aSO->GetFatherComponent();
 381       CORBA::String_var myFatherCpntDataType = myFatherCpnt->ComponentDataType();
 382       if (  strcmp(myFatherCpntDataType,"GEOM")==0) {
 383         MESSAGE("aSO father component is GEOM");
 384         if (!aSO->FindAttribute(anAttr, "AttributeIOR")) return S;
 385         myAttribute=SALOMEDS::AttributeIOR::_narrow(anAttr);
 386         myAttrValue=myAttribute->Value();
 387         MESSAGE("aSO IOR: "<< myAttrValue);
 388         myCorbaObj=smeshGen_i->GetORB()->string_to_object(myAttrValue);
 389         aShape = GEOM::GEOM_Object::_narrow(myCorbaObj);
 390       }
 391       if ( !aShape->_is_nil() )
 392         S=smeshGen_i->GeomObjectToShape( aShape.in() );
 393     }
 394     return S;
 395 }
 396 /////////////////////////////////////////////////////////
 397
 398 void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp, blsurf_session_t *bls)
 399 {
 400   int    _topology      = BLSURFPlugin_Hypothesis::GetDefaultTopology();
 401   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
 402   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
 403   int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
 404   double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
 405   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
 406   double _gradation     = BLSURFPlugin_Hypothesis::GetDefaultGradation();
 407   bool   _quadAllowed   = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
 408   bool   _decimesh      = BLSURFPlugin_Hypothesis::GetDefaultDecimesh();
 409   int    _verb          = BLSURFPlugin_Hypothesis::GetDefaultVerbosity();
 410
 411   if (hyp) {
 412     MESSAGE("BLSURFPlugin_BLSURF::SetParameters");
 413     _topology      = (int) hyp->GetTopology();
 414     _physicalMesh  = (int) hyp->GetPhysicalMesh();
 415     _phySize       = hyp->GetPhySize();
 416     _geometricMesh = (int) hyp->GetGeometricMesh();
 417     _angleMeshS    = hyp->GetAngleMeshS();
 418     _angleMeshC    = hyp->GetAngleMeshC();
 419     _gradation     = hyp->GetGradation();
 420     _quadAllowed   = hyp->GetQuadAllowed();
 421     _decimesh      = hyp->GetDecimesh();
 422     _verb          = hyp->GetVerbosity();
 423
 424     if ( hyp->GetPhyMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 425       blsurf_set_param(bls, "hphymin", to_string(hyp->GetPhyMin()).c_str());
 426     if ( hyp->GetPhyMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 427       blsurf_set_param(bls, "hphymax", to_string(hyp->GetPhyMax()).c_str());
 428     if ( hyp->GetGeoMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 429       blsurf_set_param(bls, "hgeomin", to_string(hyp->GetGeoMin()).c_str());
 430     if ( hyp->GetGeoMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 431       blsurf_set_param(bls, "hgeomax", to_string(hyp->GetGeoMax()).c_str());
 432
 433     const BLSURFPlugin_Hypothesis::TOptionValues & opts = hyp->GetOptionValues();
 434     BLSURFPlugin_Hypothesis::TOptionValues::const_iterator opIt;
 435     for ( opIt = opts.begin(); opIt != opts.end(); ++opIt )
 436       if ( !opIt->second.empty() ) {
 437         MESSAGE("blsurf_set_param(): " << opIt->first << " = " << opIt->second);
 438         blsurf_set_param(bls, opIt->first.c_str(), opIt->second.c_str());
 439       }
 440
 441   } else {
 442     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
 443   }
 444   _smp_phy_size = _phySize;
 445   blsurf_set_param(bls, "topo_points",       _topology > 0 ? "1" : "0");
 446   blsurf_set_param(bls, "topo_curves",       _topology > 0 ? "1" : "0");
 447   blsurf_set_param(bls, "topo_project",      _topology > 0 ? "1" : "0");
 448   blsurf_set_param(bls, "clean_boundary",    _topology > 1 ? "1" : "0");
 449   blsurf_set_param(bls, "close_boundary",    _topology > 1 ? "1" : "0");
 450   blsurf_set_param(bls, "hphy_flag",         to_string(_physicalMesh).c_str());
 451 //  blsurf_set_param(bls, "hphy_flag",         "2");
 452   if ((to_string(_physicalMesh))=="2"){
 453
 454     TopoDS_Shape GeomShape;
 455     TopAbs_ShapeEnum GeomType;
 456     //
 457     // Standard Size Maps
 458     //
 459     MESSAGE("Setting a Size Map");
 460     const BLSURFPlugin_Hypothesis::TSizeMap & sizeMaps = hyp->GetSizeMapEntries();
 461     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator smIt;
 462     int i=0;
 463     for ( smIt = sizeMaps.begin(); smIt != sizeMaps.end(); ++smIt ) {
 464       if ( !smIt->second.empty() ) {
 465         MESSAGE("blsurf_set_sizeMap(): " << smIt->first << " = " << smIt->second);
 466         GeomShape = entryToShape(smIt->first);
 467         GeomType  = GeomShape.ShapeType();
 468         if (GeomType == TopAbs_FACE){
 469           HasSizeMapOnFace = true;
 470           FaceId2SizeMap[TopoDS::Face(GeomShape).HashCode(471662)] = smIt->second;
 471         }
 472         if (GeomType == TopAbs_EDGE){
 473           HasSizeMapOnEdge = true;
 474           HasSizeMapOnFace = true;
 475           EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(471662)] = smIt->second;
 476         }
 477         if (GeomType == TopAbs_VERTEX){
 478           HasSizeMapOnVertex = true;
 479           HasSizeMapOnEdge   = true;
 480           HasSizeMapOnFace   = true;
 481           VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(471662)] = smIt->second;
 482         }
 483       }
 484     }
 485
 486     //
 487     // Attractors
 488     //
 489     MESSAGE("Setting Attractors");
 490     const BLSURFPlugin_Hypothesis::TSizeMap & attractors = hyp->GetAttractorEntries();
 491     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator atIt;
 492     for ( atIt = attractors.begin(); atIt != attractors.end(); ++atIt ) {
 493       if ( !atIt->second.empty() ) {
 494         MESSAGE("blsurf_set_attractor(): " << atIt->first << " = " << atIt->second);
 495         GeomShape = entryToShape(atIt->first);
 496         GeomType  = GeomShape.ShapeType();
 497
 498         if (GeomType == TopAbs_FACE){
 499           HasSizeMapOnFace = true;
 500
 501           double xa, ya, za; // Coordinates of attractor point
 502           double a, b;       // Attractor parameter
 503           int pos1, pos2;
 504           // atIt->second has the following pattern:
 505           // ATTRACTOR(xa;ya;za;a;b)
 506           // where:
 507           // xa;ya;za : coordinates of  attractor
 508           // a        : desired size on attractor
 509           // b        : distance of influence of attractor
 510           //
 511           // We search the parameters in the string
 512           pos1 = atIt->second.find(";");
 513           xa = atof(atIt->second.substr(10, pos1-10).c_str());
 514           pos2 = atIt->second.find(";", pos1+1);
 515           ya = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 516           pos1 = pos2;
 517           pos2 = atIt->second.find(";", pos1+1);
 518           za = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 519           pos1 = pos2;
 520           pos2 = atIt->second.find(";", pos1+1);
 521           a = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 522           pos1 = pos2;
 523           pos2 = atIt->second.find(")");
 524           b = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 525
 526           // Get the (u,v) values of the attractor on the face
 527           gp_XY uvPoint = getUV(TopoDS::Face(GeomShape),gp_XYZ(xa,ya,za));
 528           Standard_Real u0 = uvPoint.X();
 529           Standard_Real v0 = uvPoint.Y();
 530           // We construct the python function
 531           ostringstream attractorFunction;
 532           attractorFunction << "def f(u,v): return ";
 533           attractorFunction << _smp_phy_size << "-(" << _smp_phy_size <<"-" << a << ")";
 534           attractorFunction << "*exp(-((u-("<<u0<<"))*(u-("<<u0<<"))+(v-("<<v0<<"))*(v-("<<v0<<")))/(" << b << "*" << b <<"))";
 535
 536           MESSAGE("Python function for attractor:" << std::endl << attractorFunction.str());
 537
 538           FaceId2SizeMap[TopoDS::Face(GeomShape).HashCode(471662)] =attractorFunction.str();
 539         }
 540 /*
 541         if (GeomType == TopAbs_EDGE){
 542           HasSizeMapOnEdge = true;
 543           HasSizeMapOnFace = true;
 544         EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(471662)] = atIt->second;
 545         }
 546         if (GeomType == TopAbs_VERTEX){
 547           HasSizeMapOnVertex = true;
 548           HasSizeMapOnEdge   = true;
 549           HasSizeMapOnFace   = true;
 550         VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(471662)] = atIt->second;
 551         }
 552 */
 553       }
 554     }
 555
 556
 557 //    if (HasSizeMapOnFace){
 558     // In all size map cases (hphy_flag = 2), at least map on face must be defined
 559     MESSAGE("Setting Size Map on FACES ");
 560     blsurf_data_set_sizemap_iso_cad_face(bls, size_on_surface, &_smp_phy_size);
 561 //    }
 562
 563     if (HasSizeMapOnEdge){
 564       MESSAGE("Setting Size Map on EDGES ");
 565       blsurf_data_set_sizemap_iso_cad_edge(bls, size_on_edge, &_smp_phy_size);
 566     }
 567     if (HasSizeMapOnVertex){
 568       MESSAGE("Setting Size Map on VERTICES ");
 569       blsurf_data_set_sizemap_iso_cad_point(bls, size_on_vertex, &_smp_phy_size);
 570     }
 571   }
 572   blsurf_set_param(bls, "hphydef",           to_string(_phySize).c_str());
 573   blsurf_set_param(bls, "hgeo_flag",         to_string(_geometricMesh).c_str());
 574   blsurf_set_param(bls, "angle_meshs",       to_string(_angleMeshS).c_str());
 575   blsurf_set_param(bls, "angle_meshc",       to_string(_angleMeshC).c_str());
 576   blsurf_set_param(bls, "gradation",         to_string(_gradation).c_str());
 577   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 578   blsurf_set_param(bls, "element",           _quadAllowed ? "q1.0" : "p1");
 579   blsurf_set_param(bls, "verb",              to_string(_verb).c_str());
 580 }
 581
 582 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 583 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 584                   real *duu, real *duv, real *dvv, void *user_data);
 585 status_t message_callback(message_t *msg, void *user_data);
 586
 587 //=============================================================================
 588 /*!
 589  *
 590  */
 591 //=============================================================================
 592
 593 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 594
 595   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 596
 597   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 598     MESSAGE("  the shape is a COMPOUND");
 599   }
 600   else {
 601     MESSAGE("  the shape is UNKNOWN");
 602   };
 603
 604   context_t *ctx =  context_new();
 605   context_set_message_callback(ctx, message_callback, &_comment);
 606
 607   cad_t *c = cad_new(ctx);
 608
 609   blsurf_session_t *bls = blsurf_session_new(ctx);
 610
 611
 612   SetParameters(_hypothesis, bls);
 613
 614   TopTools_IndexedMapOfShape fmap;
 615   TopTools_IndexedMapOfShape emap;
 616   TopTools_IndexedMapOfShape pmap;
 617   vector<Handle(Geom2d_Curve)> curves;
 618   vector<Handle(Geom_Surface)> surfaces;
 619
 620
 621
 622   fmap.Clear();
 623   FaceId2PythonSmp.clear();
 624   emap.Clear();
 625   EdgeId2PythonSmp.clear();
 626   pmap.Clear();
 627   VertexId2PythonSmp.clear();
 628   surfaces.resize(0);
 629   curves.resize(0);
 630
 631   assert(Py_IsInitialized());
 632   PyGILState_STATE gstate;
 633   gstate = PyGILState_Ensure();
 634 /*
 635   Standard_Real u_min;
 636   Standard_Real v_min;
 637   Standard_Real u_max;
 638   Standard_Real v_max;
 639 */
 640   int iface = 0;
 641   string bad_end = "return";
 642   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 643     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 644     if (fmap.FindIndex(f) > 0)
 645       continue;
 646
 647     fmap.Add(f);
 648     iface++;
 649     surfaces.push_back(BRep_Tool::Surface(f));
 650     // Get bound values of uv surface
 651     //BRep_Tool::Surface(f)->Bounds(u_min,u_max,v_min,v_max);
 652     //MESSAGE("BRep_Tool::Surface(f)->Bounds(u_min,u_max,v_min,v_max): " << u_min << ", " << u_max << ", " << v_min << ", " << v_max);
 653
 654     if ((HasSizeMapOnFace) && FaceId2SizeMap.find(f.HashCode(471662))!=FaceId2SizeMap.end()){
 655         MESSAGE("FaceId2SizeMap[f.HashCode(471662)].find(bad_end): " << FaceId2SizeMap[f.HashCode(471662)].find(bad_end));
 656         MESSAGE("FaceId2SizeMap[f.HashCode(471662)].size(): " << FaceId2SizeMap[f.HashCode(471662)].size());
 657         MESSAGE("bad_end.size(): " << bad_end.size());
 658       // check if function ends with "return"
 659         if (FaceId2SizeMap[f.HashCode(471662)].find(bad_end) == (FaceId2SizeMap[f.HashCode(471662)].size()-bad_end.size()-1))
 660         continue;
 661       // Expr To Python function, verification is performed at validation in GUI
 662       PyObject * obj = NULL;
 663       obj= PyRun_String(FaceId2SizeMap[f.HashCode(471662)].c_str(), Py_file_input, main_dict, NULL);
 664       Py_DECREF(obj);
 665       PyObject * func = NULL;
 666       func = PyObject_GetAttrString(main_mod, "f");
 667       FaceId2PythonSmp[iface]=func;
 668       FaceId2SizeMap.erase(f.HashCode(471662));
 669     }
 670     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 671     cad_face_set_tag(fce, iface);
 672     if(f.Orientation() != TopAbs_FORWARD){
 673       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 674     } else {
 675       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 676     }
 677
 678     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
 679       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
 680       int ic = emap.FindIndex(e);
 681       if (ic <= 0)
 682         ic = emap.Add(e);
 683
 684       double tmin,tmax;
 685       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
 686       if ((HasSizeMapOnEdge) && EdgeId2SizeMap.find(e.HashCode(471662))!=EdgeId2SizeMap.end()){
 687           if (EdgeId2SizeMap[e.HashCode(471662)].find(bad_end) == (EdgeId2SizeMap[e.HashCode(471662)].size()-bad_end.size()-1))
 688           continue;
 689         // Expr To Python function, verification is performed at validation in GUI
 690         PyObject * obj = NULL;
 691         obj= PyRun_String(EdgeId2SizeMap[e.HashCode(471662)].c_str(), Py_file_input, main_dict, NULL);
 692         Py_DECREF(obj);
 693         PyObject * func = NULL;
 694         func = PyObject_GetAttrString(main_mod, "f");
 695         EdgeId2PythonSmp[ic]=func;
 696         EdgeId2SizeMap.erase(e.HashCode(471662));
 697       }
 698       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
 699       cad_edge_set_tag(edg, ic);
 700       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
 701       if (e.Orientation() == TopAbs_INTERNAL)
 702         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
 703
 704       int npts = 0;
 705       int ip1, ip2, *ip;
 706       gp_Pnt2d e0 = curves.back()->Value(tmin);
 707       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
 708       Standard_Real d1=0,d2=0;
 709       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
 710         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
 711         ++npts;
 712         if (npts == 1){
 713           ip = &ip1;
 714           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 715         } else {
 716           ip = &ip2;
 717           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 718         }
 719         *ip = pmap.FindIndex(v);
 720         if(*ip <= 0)
 721           *ip = pmap.Add(v);
 722     if ((HasSizeMapOnVertex) && VertexId2SizeMap.find(v.HashCode(471662))!=VertexId2SizeMap.end()){
 723         if (VertexId2SizeMap[v.HashCode(471662)].find(bad_end) == (VertexId2SizeMap[v.HashCode(471662)].size()-bad_end.size()-1))
 724             continue;
 725           // Expr To Python function, verification is performed at validation in GUI
 726           PyObject * obj = NULL;
 727           obj= PyRun_String(VertexId2SizeMap[v.HashCode(471662)].c_str(), Py_file_input, main_dict, NULL);
 728           Py_DECREF(obj);
 729           PyObject * func = NULL;
 730           func = PyObject_GetAttrString(main_mod, "f");
 731           VertexId2PythonSmp[*ip]=func;
 732           VertexId2SizeMap.erase(v.HashCode(471662));
 733         }
 734       }
 735       if (npts != 2) {
 736         // should not happen
 737         MESSAGE("An edge does not have 2 extremities.");
 738       } else {
 739         if (d1 < d2)
 740           cad_edge_set_extremities(edg, ip1, ip2);
 741         else
 742           cad_edge_set_extremities(edg, ip2, ip1);
 743       }
 744     } // for edge
 745   } //for face
 746
 747
 748   PyGILState_Release(gstate);
 749
 750   blsurf_data_set_cad(bls, c);
 751
 752   std::cout << std::endl;
 753   std::cout << "Beginning of Surface Mesh generation" << std::endl;
 754   std::cout << std::endl;
 755
 756   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 757 #ifndef WNT
 758   feclearexcept( FE_ALL_EXCEPT );
 759   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
 760 #endif
 761
 762     status_t status = STATUS_ERROR;
 763
 764   try {
 765     OCC_CATCH_SIGNALS;
 766
 767     status = blsurf_compute_mesh(bls);
 768
 769   }
 770   catch ( std::exception& exc ) {
 771     _comment += exc.what();
 772   }
 773   catch (Standard_Failure& ex) {
 774     _comment += ex.DynamicType()->Name();
 775     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
 776       _comment += ": ";
 777       _comment += ex.GetMessageString();
 778     }
 779   }
 780   catch (...) {
 781     if ( _comment.empty() )
 782       _comment = "Exception in blsurf_compute_mesh()";
 783   }
 784   if ( status != STATUS_OK) {
 785     blsurf_session_delete(bls);
 786     cad_delete(c);
 787     context_delete(ctx);
 788
 789     return error(_comment);
 790   }
 791
 792   std::cout << std::endl;
 793   std::cout << "End of Surface Mesh generation" << std::endl;
 794   std::cout << std::endl;
 795
 796   mesh_t *msh;
 797   blsurf_data_get_mesh(bls, &msh);
 798   if(!msh){
 799     blsurf_session_delete(bls);
 800     cad_delete(c);
 801     context_delete(ctx);
 802
 803     return error(_comment);
 804     //return false;
 805   }
 806
 807   integer nv, ne, nt, nq, vtx[4], tag;
 808   real xyz[3];
 809
 810   mesh_get_vertex_count(msh, &nv);
 811   mesh_get_edge_count(msh, &ne);
 812   mesh_get_triangle_count(msh, &nt);
 813   mesh_get_quadrangle_count(msh, &nq);
 814
 815
 816   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 817   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
 818   bool* tags = new bool[nv+1];
 819
 820   for(int iv=1;iv<=nv;iv++) {
 821     mesh_get_vertex_coordinates(msh, iv, xyz);
 822     mesh_get_vertex_tag(msh, iv, &tag);
 823     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
 824     // internal point are tagged to zero
 825     if(tag){
 826       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
 827       tags[iv] = false;
 828     } else {
 829       tags[iv] = true;
 830     }
 831   }
 832
 833   for(int it=1;it<=ne;it++) {
 834     mesh_get_edge_vertices(msh, it, vtx);
 835     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
 836     mesh_get_edge_tag(msh, it, &tag);
 837
 838     if (tags[vtx[0]]) {
 839       meshDS->SetNodeOnEdge(nodes[vtx[0]], TopoDS::Edge(emap(tag)));
 840       tags[vtx[0]] = false;
 841     };
 842     if (tags[vtx[1]]) {
 843       meshDS->SetNodeOnEdge(nodes[vtx[1]], TopoDS::Edge(emap(tag)));
 844       tags[vtx[1]] = false;
 845     };
 846     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
 847
 848   }
 849
 850   for(int it=1;it<=nt;it++) {
 851     mesh_get_triangle_vertices(msh, it, vtx);
 852     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
 853     mesh_get_triangle_tag(msh, it, &tag);
 854     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
 855     if (tags[vtx[0]]) {
 856       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 857       tags[vtx[0]] = false;
 858     };
 859     if (tags[vtx[1]]) {
 860       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 861       tags[vtx[1]] = false;
 862     };
 863     if (tags[vtx[2]]) {
 864       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 865       tags[vtx[2]] = false;
 866     };
 867   }
 868
 869   for(int it=1;it<=nq;it++) {
 870     mesh_get_quadrangle_vertices(msh, it, vtx);
 871     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
 872     mesh_get_quadrangle_tag(msh, it, &tag);
 873     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
 874     if (tags[vtx[0]]) {
 875       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 876       tags[vtx[0]] = false;
 877     };
 878     if (tags[vtx[1]]) {
 879       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 880       tags[vtx[1]] = false;
 881     };
 882     if (tags[vtx[2]]) {
 883       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 884       tags[vtx[2]] = false;
 885     };
 886     if (tags[vtx[3]]) {
 887       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
 888       tags[vtx[3]] = false;
 889     };
 890   }
 891
 892   delete nodes;
 893
 894   /* release the mesh object */
 895   blsurf_data_regain_mesh(bls, msh);
 896
 897   /* clean up everything */
 898   blsurf_session_delete(bls);
 899   cad_delete(c);
 900
 901   context_delete(ctx);
 902
 903   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 904 #ifndef WNT
 905   if ( oldFEFlags > 0 )
 906     feenableexcept( oldFEFlags );
 907   feclearexcept( FE_ALL_EXCEPT );
 908 #endif
 909
 910   return true;
 911 }
 912
 913 //=============================================================================
 914 /*!
 915  *
 916  */
 917 //=============================================================================
 918
 919 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
 920 {
 921   return save;
 922 }
 923
 924 //=============================================================================
 925 /*!
 926  *
 927  */
 928 //=============================================================================
 929
 930 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
 931 {
 932   return load;
 933 }
 934
 935 //=============================================================================
 936 /*!
 937  *
 938  */
 939 //=============================================================================
 940
 941 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
 942 {
 943   return hyp.SaveTo( save );
 944 }
 945
 946 //=============================================================================
 947 /*!
 948  *
 949  */
 950 //=============================================================================
 951
 952 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
 953 {
 954   return hyp.LoadFrom( load );
 955 }
 956
 957 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
 958 {
 959   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
 960
 961   if (uv){
 962     gp_Pnt2d P;
 963     P=pargeo->Value(t);
 964     uv[0]=P.X(); uv[1]=P.Y();
 965   }
 966
 967   if(dt) {
 968     gp_Vec2d V1;
 969     V1=pargeo->DN(t,1);
 970     dt[0]=V1.X(); dt[1]=V1.Y();
 971   }
 972
 973   if(dtt){
 974     gp_Vec2d V2;
 975     V2=pargeo->DN(t,2);
 976     dtt[0]=V2.X(); dtt[1]=V2.Y();
 977   }
 978
 979   return 0;
 980 }
 981
 982 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 983                   real *duu, real *duv, real *dvv, void *user_data)
 984 {
 985   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
 986
 987   if(xyz){
 988    gp_Pnt P;
 989    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
 990    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
 991   }
 992
 993   if(du && dv){
 994     gp_Pnt P;
 995     gp_Vec D1U,D1V;
 996
 997     geometry->D1(uv[0],uv[1],P,D1U,D1V);
 998     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
 999     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
1000   }
1001
1002   if(duu && duv && dvv){
1003     gp_Pnt P;
1004     gp_Vec D1U,D1V;
1005     gp_Vec D2U,D2V,D2UV;
1006
1007     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
1008     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
1009     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
1010     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
1011   }
1012
1013   return 0;
1014 }
1015
1016
1017 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
1018 {
1019   if (face_id == 1) {
1020     if (my_u_min > uv[0]) {
1021       my_u_min = uv[0];
1022     }
1023     if (my_v_min > uv[1]) {
1024       my_v_min = uv[1];
1025     }
1026     if (my_u_max < uv[0]) {
1027       my_u_max = uv[0];
1028     }
1029     if (my_v_max < uv[1]) {
1030       my_v_max = uv[1];
1031     }
1032   }
1033
1034   if (FaceId2PythonSmp.count(face_id) != 0){
1035     PyObject * pyresult = NULL;
1036     PyObject* new_stderr = NULL;
1037     assert(Py_IsInitialized());
1038     PyGILState_STATE gstate;
1039     gstate = PyGILState_Ensure();
1040     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],"(f,f)",uv[0],uv[1]);
1041     double result;
1042     if ( pyresult == NULL){
1043       fflush(stderr);
1044       string err_description="";
1045       new_stderr = newPyStdOut(err_description);
1046       PySys_SetObject("stderr", new_stderr);
1047       PyErr_Print();
1048       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1049       Py_DECREF(new_stderr);
1050       MESSAGE("Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description);
1051       result = *((double*)user_data);
1052       }
1053     else {
1054       result = PyFloat_AsDouble(pyresult);
1055       Py_DECREF(pyresult);
1056     }
1057     *size = result;
1058     //MESSAGE("f(" << uv[0] << "," << uv[1] << ")" << " = " << result);
1059     PyGILState_Release(gstate);
1060   }
1061   else {
1062     *size = *((double*)user_data);
1063   }
1064   return STATUS_OK;
1065 }
1066
1067 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
1068 {
1069   if (EdgeId2PythonSmp.count(edge_id) != 0){
1070     PyObject * pyresult = NULL;
1071     PyObject* new_stderr = NULL;
1072     assert(Py_IsInitialized());
1073     PyGILState_STATE gstate;
1074     gstate = PyGILState_Ensure();
1075     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],"(f)",t);
1076     double result;
1077     if ( pyresult == NULL){
1078       fflush(stderr);
1079       string err_description="";
1080       new_stderr = newPyStdOut(err_description);
1081       PySys_SetObject("stderr", new_stderr);
1082       PyErr_Print();
1083       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1084       Py_DECREF(new_stderr);
1085       MESSAGE("Can't evaluate f(" << t << ")" << " error is " << err_description);
1086       result = *((double*)user_data);
1087       }
1088     else {
1089       result = PyFloat_AsDouble(pyresult);
1090       Py_DECREF(pyresult);
1091     }
1092     *size = result;
1093     PyGILState_Release(gstate);
1094   }
1095   else {
1096     *size = *((double*)user_data);
1097   }
1098   return STATUS_OK;
1099 }
1100
1101 status_t size_on_vertex(integer point_id, real *size, void *user_data)
1102 {
1103   if (VertexId2PythonSmp.count(point_id) != 0){
1104     PyObject * pyresult = NULL;
1105     PyObject* new_stderr = NULL;
1106     assert(Py_IsInitialized());
1107     PyGILState_STATE gstate;
1108     gstate = PyGILState_Ensure();
1109     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],"");
1110     double result;
1111     if ( pyresult == NULL){
1112       fflush(stderr);
1113       string err_description="";
1114       new_stderr = newPyStdOut(err_description);
1115       PySys_SetObject("stderr", new_stderr);
1116       PyErr_Print();
1117       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1118       Py_DECREF(new_stderr);
1119       MESSAGE("Can't evaluate f()" << " error is " << err_description);
1120       result = *((double*)user_data);
1121       }
1122     else {
1123       result = PyFloat_AsDouble(pyresult);
1124       Py_DECREF(pyresult);
1125     }
1126     *size = result;
1127     PyGILState_Release(gstate);
1128   }
1129   else {
1130     *size = *((double*)user_data);
1131   }
1132  return STATUS_OK;
1133 }
1134
1135 status_t message_callback(message_t *msg, void *user_data)
1136 {
1137   integer errnumber = 0;
1138   char *desc;
1139   message_get_number(msg, &errnumber);
1140   message_get_description(msg, &desc);
1141   if ( errnumber < 0 ) {
1142     string * error = (string*)user_data;
1143 //   if ( !error->empty() )
1144 //     *error += "\n";
1145     // remove ^A from the tail
1146     int len = strlen( desc );
1147     while (len > 0 && desc[len-1] != '\n')
1148       len--;
1149     error->append( desc, len );
1150   }
1151   else {
1152       std::cout << desc << std::endl;
1153   }
1154   return STATUS_OK;
1155 }
1156
1157
1158 //=============================================================================
1159 /*!
1160  *
1161  */
1162 //=============================================================================
1163 bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh& aMesh,
1164                                    const TopoDS_Shape& aShape,
1165                                    MapShapeNbElems& aResMap)
1166 {
1167   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
1168   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
1169   //int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
1170   //double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
1171   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
1172   if(_hypothesis) {
1173     _physicalMesh  = (int) _hypothesis->GetPhysicalMesh();
1174     _phySize       = _hypothesis->GetPhySize();
1175     //_geometricMesh = (int) hyp->GetGeometricMesh();
1176     //_angleMeshS    = hyp->GetAngleMeshS();
1177     _angleMeshC    = _hypothesis->GetAngleMeshC();
1178   }
1179
1180   bool IsQuadratic = false;
1181
1182   // ----------------
1183   // evaluate 1D
1184   // ----------------
1185   TopTools_DataMapOfShapeInteger EdgesMap;
1186   double fullLen = 0.0;
1187   double fullNbSeg = 0;
1188   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
1189     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1190     if( EdgesMap.IsBound(E) )
1191       continue;
1192     SMESH_subMesh *sm = aMesh.GetSubMesh(E);
1193     double aLen = SMESH_Algo::EdgeLength(E);
1194     fullLen += aLen;
1195     int nb1d = 0;
1196     if(_physicalMesh==1) {
1197        nb1d = (int)aLen/_phySize + 1;
1198     }
1199     else {
1200       // use geometry
1201       double f,l;
1202       Handle(Geom_Curve) C = BRep_Tool::Curve(E,f,l);
1203       double fullAng = 0.0;
1204       double dp = (l-f)/200;
1205       gp_Pnt P1,P2,P3;
1206       C->D0(f,P1);
1207       C->D0(f+dp,P2);
1208       gp_Vec V1(P1,P2);
1209       for(int j=2; j<=200; j++) {
1210         C->D0(f+dp*j,P3);
1211         gp_Vec V2(P2,P3);
1212         fullAng += fabs(V1.Angle(V2));
1213         V1 = V2;
1214         P2 = P3;
1215       }
1216       nb1d = fullAng/_angleMeshC + 1;
1217     }
1218     fullNbSeg += nb1d;
1219     std::vector<int> aVec(17);
1220     for(int i=0; i<17; i++) aVec[i]=0;
1221     if( IsQuadratic > 0 ) {
1222       aVec[0] = 2*nb1d - 1;
1223       aVec[2] = nb1d;
1224     }
1225     else {
1226       aVec[0] = nb1d - 1;
1227       aVec[1] = nb1d;
1228     }
1229     aResMap.insert(std::make_pair(sm,aVec));
1230     EdgesMap.Bind(E,nb1d);
1231   }
1232   double ELen = fullLen/fullNbSeg;
1233   // ----------------
1234   // evaluate 2D
1235   // ----------------
1236   // try to evaluate as in MEFISTO
1237   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
1238     TopoDS_Face F = TopoDS::Face( exp.Current() );
1239     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1240     GProp_GProps G;
1241     BRepGProp::SurfaceProperties(F,G);
1242     double anArea = G.Mass();
1243     int nb1d = 0;
1244     for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) {
1245       nb1d += EdgesMap.Find(exp1.Current());
1246     }
1247     int nbFaces = (int) anArea/(ELen*ELen*sqrt(3)/4);
1248     int nbNodes = (int) ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1;
1249     std::vector<int> aVec(17);
1250     for(int i=0; i<17; i++) aVec[i]=0;
1251     if( IsQuadratic ) {
1252       int nb1d_in = (nbFaces*3 - nb1d) / 2;
1253       aVec[0] = nbNodes + nb1d_in;
1254       aVec[4] = nbFaces;
1255     }
1256     else {
1257       aVec[0] = nbNodes;
1258       aVec[3] = nbFaces;
1259     }
1260     aResMap.insert(std::make_pair(sm,aVec));
1261   }
1262
1263   // ----------------
1264   // evaluate 3D
1265   // ----------------
1266   GProp_GProps G;
1267   BRepGProp::VolumeProperties(aShape,G);
1268   double aVolume = G.Mass();
1269   double tetrVol = 0.1179*ELen*ELen*ELen;
1270   int nbVols = (int)aVolume/tetrVol;
1271   int nb1d_in = (int) ( nbVols*6 - fullNbSeg ) / 6;
1272   std::vector<int> aVec(17);
1273   for(int i=0; i<17; i++) aVec[i]=0;
1274   if( IsQuadratic ) {
1275     aVec[0] = nb1d_in/3 + 1 + nb1d_in;
1276     aVec[9] = nbVols;
1277   }
1278   else {
1279     aVec[0] = nb1d_in/3 + 1;
1280     aVec[8] = nbVols;
1281   }
1282   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1283   aResMap.insert(std::make_pair(sm,aVec));
1284
1285   return true;
1286 }