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