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     for ( smIt = sizeMaps.begin(); smIt != sizeMaps.end(); ++smIt ) {
 463       if ( !smIt->second.empty() ) {
 464         MESSAGE("blsurf_set_sizeMap(): " << smIt->first << " = " << smIt->second);
 465         GeomShape = entryToShape(smIt->first);
 466         GeomType  = GeomShape.ShapeType();
 467         if (GeomType == TopAbs_FACE){
 468           HasSizeMapOnFace = true;
 469           FaceId2SizeMap[TopoDS::Face(GeomShape).HashCode(471662)] = smIt->second;
 470         }
 471         if (GeomType == TopAbs_EDGE){
 472           HasSizeMapOnEdge = true;
 473           HasSizeMapOnFace = true;
 474           EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(471662)] = smIt->second;
 475         }
 476         if (GeomType == TopAbs_VERTEX){
 477           HasSizeMapOnVertex = true;
 478           HasSizeMapOnEdge   = true;
 479           HasSizeMapOnFace   = true;
 480           VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(471662)] = smIt->second;
 481         }
 482       }
 483     }
 484
 485     //
 486     // Attractors
 487     //
 488     MESSAGE("Setting Attractors");
 489     const BLSURFPlugin_Hypothesis::TSizeMap & attractors = hyp->GetAttractorEntries();
 490     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator atIt;
 491     for ( atIt = attractors.begin(); atIt != attractors.end(); ++atIt ) {
 492       if ( !atIt->second.empty() ) {
 493         MESSAGE("blsurf_set_attractor(): " << atIt->first << " = " << atIt->second);
 494         GeomShape = entryToShape(atIt->first);
 495         GeomType  = GeomShape.ShapeType();
 496
 497         if (GeomType == TopAbs_FACE){
 498           HasSizeMapOnFace = true;
 499
 500           double xa, ya, za; // Coordinates of attractor point
 501           double a, b;       // Attractor parameter
 502           int pos1, pos2;
 503           // atIt->second has the following pattern:
 504           // ATTRACTOR(xa;ya;za;a;b)
 505           // where:
 506           // xa;ya;za : coordinates of  attractor
 507           // a        : desired size on attractor
 508           // b        : distance of influence of attractor
 509           //
 510           // We search the parameters in the string
 511           pos1 = atIt->second.find(";");
 512           xa = atof(atIt->second.substr(10, pos1-10).c_str());
 513           pos2 = atIt->second.find(";", pos1+1);
 514           ya = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 515           pos1 = pos2;
 516           pos2 = atIt->second.find(";", pos1+1);
 517           za = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 518           pos1 = pos2;
 519           pos2 = atIt->second.find(";", pos1+1);
 520           a = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 521           pos1 = pos2;
 522           pos2 = atIt->second.find(")");
 523           b = atof(atIt->second.substr(pos1+1, pos2-pos1-1).c_str());
 524
 525           // Get the (u,v) values of the attractor on the face
 526           gp_XY uvPoint = getUV(TopoDS::Face(GeomShape),gp_XYZ(xa,ya,za));
 527           Standard_Real u0 = uvPoint.X();
 528           Standard_Real v0 = uvPoint.Y();
 529           // We construct the python function
 530           ostringstream attractorFunction;
 531           attractorFunction << "def f(u,v): return ";
 532           attractorFunction << _smp_phy_size << "-(" << _smp_phy_size <<"-" << a << ")";
 533           attractorFunction << "*exp(-((u-("<<u0<<"))*(u-("<<u0<<"))+(v-("<<v0<<"))*(v-("<<v0<<")))/(" << b << "*" << b <<"))";
 534
 535           MESSAGE("Python function for attractor:" << std::endl << attractorFunction.str());
 536
 537           FaceId2SizeMap[TopoDS::Face(GeomShape).HashCode(471662)] =attractorFunction.str();
 538         }
 539 /*
 540         if (GeomType == TopAbs_EDGE){
 541           HasSizeMapOnEdge = true;
 542           HasSizeMapOnFace = true;
 543         EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(471662)] = atIt->second;
 544         }
 545         if (GeomType == TopAbs_VERTEX){
 546           HasSizeMapOnVertex = true;
 547           HasSizeMapOnEdge   = true;
 548           HasSizeMapOnFace   = true;
 549         VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(471662)] = atIt->second;
 550         }
 551 */
 552       }
 553     }
 554
 555
 556 //    if (HasSizeMapOnFace){
 557     // In all size map cases (hphy_flag = 2), at least map on face must be defined
 558     MESSAGE("Setting Size Map on FACES ");
 559     blsurf_data_set_sizemap_iso_cad_face(bls, size_on_surface, &_smp_phy_size);
 560 //    }
 561
 562     if (HasSizeMapOnEdge){
 563       MESSAGE("Setting Size Map on EDGES ");
 564       blsurf_data_set_sizemap_iso_cad_edge(bls, size_on_edge, &_smp_phy_size);
 565     }
 566     if (HasSizeMapOnVertex){
 567       MESSAGE("Setting Size Map on VERTICES ");
 568       blsurf_data_set_sizemap_iso_cad_point(bls, size_on_vertex, &_smp_phy_size);
 569     }
 570   }
 571   blsurf_set_param(bls, "hphydef",           to_string(_phySize).c_str());
 572   blsurf_set_param(bls, "hgeo_flag",         to_string(_geometricMesh).c_str());
 573   blsurf_set_param(bls, "angle_meshs",       to_string(_angleMeshS).c_str());
 574   blsurf_set_param(bls, "angle_meshc",       to_string(_angleMeshC).c_str());
 575   blsurf_set_param(bls, "gradation",         to_string(_gradation).c_str());
 576   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 577   blsurf_set_param(bls, "element",           _quadAllowed ? "q1.0" : "p1");
 578   blsurf_set_param(bls, "verb",              to_string(_verb).c_str());
 579 }
 580
 581 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 582 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 583                   real *duu, real *duv, real *dvv, void *user_data);
 584 status_t message_callback(message_t *msg, void *user_data);
 585
 586 //=============================================================================
 587 /*!
 588  *
 589  */
 590 //=============================================================================
 591
 592 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 593
 594   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 595
 596   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 597     MESSAGE("  the shape is a COMPOUND");
 598   }
 599   else {
 600     MESSAGE("  the shape is UNKNOWN");
 601   };
 602
 603   context_t *ctx =  context_new();
 604   context_set_message_callback(ctx, message_callback, &_comment);
 605
 606   cad_t *c = cad_new(ctx);
 607
 608   blsurf_session_t *bls = blsurf_session_new(ctx);
 609
 610
 611   SetParameters(_hypothesis, bls);
 612
 613   TopTools_IndexedMapOfShape fmap;
 614   TopTools_IndexedMapOfShape emap;
 615   TopTools_IndexedMapOfShape pmap;
 616   vector<Handle(Geom2d_Curve)> curves;
 617   vector<Handle(Geom_Surface)> surfaces;
 618
 619
 620
 621   fmap.Clear();
 622   FaceId2PythonSmp.clear();
 623   emap.Clear();
 624   EdgeId2PythonSmp.clear();
 625   pmap.Clear();
 626   VertexId2PythonSmp.clear();
 627   surfaces.resize(0);
 628   curves.resize(0);
 629
 630   assert(Py_IsInitialized());
 631   PyGILState_STATE gstate;
 632   gstate = PyGILState_Ensure();
 633 /*
 634   Standard_Real u_min;
 635   Standard_Real v_min;
 636   Standard_Real u_max;
 637   Standard_Real v_max;
 638 */
 639   int iface = 0;
 640   string bad_end = "return";
 641   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 642     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 643     if (fmap.FindIndex(f) > 0)
 644       continue;
 645
 646     fmap.Add(f);
 647     iface++;
 648     surfaces.push_back(BRep_Tool::Surface(f));
 649     // Get bound values of uv surface
 650     //BRep_Tool::Surface(f)->Bounds(u_min,u_max,v_min,v_max);
 651     //MESSAGE("BRep_Tool::Surface(f)->Bounds(u_min,u_max,v_min,v_max): " << u_min << ", " << u_max << ", " << v_min << ", " << v_max);
 652
 653     if ((HasSizeMapOnFace) && FaceId2SizeMap.find(f.HashCode(471662))!=FaceId2SizeMap.end()){
 654         MESSAGE("FaceId2SizeMap[f.HashCode(471662)].find(bad_end): " << FaceId2SizeMap[f.HashCode(471662)].find(bad_end));
 655         MESSAGE("FaceId2SizeMap[f.HashCode(471662)].size(): " << FaceId2SizeMap[f.HashCode(471662)].size());
 656         MESSAGE("bad_end.size(): " << bad_end.size());
 657       // check if function ends with "return"
 658         if (FaceId2SizeMap[f.HashCode(471662)].find(bad_end) == (FaceId2SizeMap[f.HashCode(471662)].size()-bad_end.size()-1))
 659         continue;
 660       // Expr To Python function, verification is performed at validation in GUI
 661       PyObject * obj = NULL;
 662       obj= PyRun_String(FaceId2SizeMap[f.HashCode(471662)].c_str(), Py_file_input, main_dict, NULL);
 663       Py_DECREF(obj);
 664       PyObject * func = NULL;
 665       func = PyObject_GetAttrString(main_mod, "f");
 666       FaceId2PythonSmp[iface]=func;
 667       FaceId2SizeMap.erase(f.HashCode(471662));
 668     }
 669     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 670     cad_face_set_tag(fce, iface);
 671     if(f.Orientation() != TopAbs_FORWARD){
 672       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 673     } else {
 674       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 675     }
 676
 677     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
 678       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
 679       int ic = emap.FindIndex(e);
 680       if (ic <= 0)
 681         ic = emap.Add(e);
 682
 683       double tmin,tmax;
 684       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
 685       if ((HasSizeMapOnEdge) && EdgeId2SizeMap.find(e.HashCode(471662))!=EdgeId2SizeMap.end()){
 686           if (EdgeId2SizeMap[e.HashCode(471662)].find(bad_end) == (EdgeId2SizeMap[e.HashCode(471662)].size()-bad_end.size()-1))
 687           continue;
 688         // Expr To Python function, verification is performed at validation in GUI
 689         PyObject * obj = NULL;
 690         obj= PyRun_String(EdgeId2SizeMap[e.HashCode(471662)].c_str(), Py_file_input, main_dict, NULL);
 691         Py_DECREF(obj);
 692         PyObject * func = NULL;
 693         func = PyObject_GetAttrString(main_mod, "f");
 694         EdgeId2PythonSmp[ic]=func;
 695         EdgeId2SizeMap.erase(e.HashCode(471662));
 696       }
 697       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
 698       cad_edge_set_tag(edg, ic);
 699       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
 700       if (e.Orientation() == TopAbs_INTERNAL)
 701         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
 702
 703       int npts = 0;
 704       int ip1, ip2, *ip;
 705       gp_Pnt2d e0 = curves.back()->Value(tmin);
 706       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
 707       Standard_Real d1=0,d2=0;
 708       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
 709         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
 710         ++npts;
 711         if (npts == 1){
 712           ip = &ip1;
 713           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 714         } else {
 715           ip = &ip2;
 716           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 717         }
 718         *ip = pmap.FindIndex(v);
 719         if(*ip <= 0)
 720           *ip = pmap.Add(v);
 721     if ((HasSizeMapOnVertex) && VertexId2SizeMap.find(v.HashCode(471662))!=VertexId2SizeMap.end()){
 722         if (VertexId2SizeMap[v.HashCode(471662)].find(bad_end) == (VertexId2SizeMap[v.HashCode(471662)].size()-bad_end.size()-1))
 723             continue;
 724           // Expr To Python function, verification is performed at validation in GUI
 725           PyObject * obj = NULL;
 726           obj= PyRun_String(VertexId2SizeMap[v.HashCode(471662)].c_str(), Py_file_input, main_dict, NULL);
 727           Py_DECREF(obj);
 728           PyObject * func = NULL;
 729           func = PyObject_GetAttrString(main_mod, "f");
 730           VertexId2PythonSmp[*ip]=func;
 731           VertexId2SizeMap.erase(v.HashCode(471662));
 732         }
 733       }
 734       if (npts != 2) {
 735         // should not happen
 736         MESSAGE("An edge does not have 2 extremities.");
 737       } else {
 738         if (d1 < d2)
 739           cad_edge_set_extremities(edg, ip1, ip2);
 740         else
 741           cad_edge_set_extremities(edg, ip2, ip1);
 742       }
 743     } // for edge
 744   } //for face
 745
 746
 747   PyGILState_Release(gstate);
 748
 749   blsurf_data_set_cad(bls, c);
 750
 751   std::cout << std::endl;
 752   std::cout << "Beginning of Surface Mesh generation" << std::endl;
 753   std::cout << std::endl;
 754
 755   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 756 #ifndef WNT
 757   feclearexcept( FE_ALL_EXCEPT );
 758   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
 759 #endif
 760
 761     status_t status = STATUS_ERROR;
 762
 763   try {
 764     OCC_CATCH_SIGNALS;
 765
 766     status = blsurf_compute_mesh(bls);
 767
 768   }
 769   catch ( std::exception& exc ) {
 770     _comment += exc.what();
 771   }
 772   catch (Standard_Failure& ex) {
 773     _comment += ex.DynamicType()->Name();
 774     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
 775       _comment += ": ";
 776       _comment += ex.GetMessageString();
 777     }
 778   }
 779   catch (...) {
 780     if ( _comment.empty() )
 781       _comment = "Exception in blsurf_compute_mesh()";
 782   }
 783   if ( status != STATUS_OK) {
 784     blsurf_session_delete(bls);
 785     cad_delete(c);
 786     context_delete(ctx);
 787
 788     return error(_comment);
 789   }
 790
 791   std::cout << std::endl;
 792   std::cout << "End of Surface Mesh generation" << std::endl;
 793   std::cout << std::endl;
 794
 795   mesh_t *msh;
 796   blsurf_data_get_mesh(bls, &msh);
 797   if(!msh){
 798     blsurf_session_delete(bls);
 799     cad_delete(c);
 800     context_delete(ctx);
 801
 802     return error(_comment);
 803     //return false;
 804   }
 805
 806   integer nv, ne, nt, nq, vtx[4], tag;
 807   real xyz[3];
 808
 809   mesh_get_vertex_count(msh, &nv);
 810   mesh_get_edge_count(msh, &ne);
 811   mesh_get_triangle_count(msh, &nt);
 812   mesh_get_quadrangle_count(msh, &nq);
 813
 814
 815   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 816   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
 817   bool* tags = new bool[nv+1];
 818
 819   for(int iv=1;iv<=nv;iv++) {
 820     mesh_get_vertex_coordinates(msh, iv, xyz);
 821     mesh_get_vertex_tag(msh, iv, &tag);
 822     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
 823     // internal point are tagged to zero
 824     if(tag){
 825       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
 826       tags[iv] = false;
 827     } else {
 828       tags[iv] = true;
 829     }
 830   }
 831
 832   for(int it=1;it<=ne;it++) {
 833     mesh_get_edge_vertices(msh, it, vtx);
 834     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
 835     mesh_get_edge_tag(msh, it, &tag);
 836
 837     if (tags[vtx[0]]) {
 838       meshDS->SetNodeOnEdge(nodes[vtx[0]], TopoDS::Edge(emap(tag)));
 839       tags[vtx[0]] = false;
 840     };
 841     if (tags[vtx[1]]) {
 842       meshDS->SetNodeOnEdge(nodes[vtx[1]], TopoDS::Edge(emap(tag)));
 843       tags[vtx[1]] = false;
 844     };
 845     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
 846
 847   }
 848
 849   for(int it=1;it<=nt;it++) {
 850     mesh_get_triangle_vertices(msh, it, vtx);
 851     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
 852     mesh_get_triangle_tag(msh, it, &tag);
 853     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
 854     if (tags[vtx[0]]) {
 855       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 856       tags[vtx[0]] = false;
 857     };
 858     if (tags[vtx[1]]) {
 859       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 860       tags[vtx[1]] = false;
 861     };
 862     if (tags[vtx[2]]) {
 863       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 864       tags[vtx[2]] = false;
 865     };
 866   }
 867
 868   for(int it=1;it<=nq;it++) {
 869     mesh_get_quadrangle_vertices(msh, it, vtx);
 870     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
 871     mesh_get_quadrangle_tag(msh, it, &tag);
 872     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
 873     if (tags[vtx[0]]) {
 874       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 875       tags[vtx[0]] = false;
 876     };
 877     if (tags[vtx[1]]) {
 878       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 879       tags[vtx[1]] = false;
 880     };
 881     if (tags[vtx[2]]) {
 882       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 883       tags[vtx[2]] = false;
 884     };
 885     if (tags[vtx[3]]) {
 886       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
 887       tags[vtx[3]] = false;
 888     };
 889   }
 890
 891   delete nodes;
 892
 893   /* release the mesh object */
 894   blsurf_data_regain_mesh(bls, msh);
 895
 896   /* clean up everything */
 897   blsurf_session_delete(bls);
 898   cad_delete(c);
 899
 900   context_delete(ctx);
 901
 902   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 903 #ifndef WNT
 904   if ( oldFEFlags > 0 )
 905     feenableexcept( oldFEFlags );
 906   feclearexcept( FE_ALL_EXCEPT );
 907 #endif
 908
 909   return true;
 910 }
 911
 912 //=============================================================================
 913 /*!
 914  *
 915  */
 916 //=============================================================================
 917
 918 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
 919 {
 920   return save;
 921 }
 922
 923 //=============================================================================
 924 /*!
 925  *
 926  */
 927 //=============================================================================
 928
 929 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
 930 {
 931   return load;
 932 }
 933
 934 //=============================================================================
 935 /*!
 936  *
 937  */
 938 //=============================================================================
 939
 940 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
 941 {
 942   return hyp.SaveTo( save );
 943 }
 944
 945 //=============================================================================
 946 /*!
 947  *
 948  */
 949 //=============================================================================
 950
 951 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
 952 {
 953   return hyp.LoadFrom( load );
 954 }
 955
 956 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
 957 {
 958   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
 959
 960   if (uv){
 961     gp_Pnt2d P;
 962     P=pargeo->Value(t);
 963     uv[0]=P.X(); uv[1]=P.Y();
 964   }
 965
 966   if(dt) {
 967     gp_Vec2d V1;
 968     V1=pargeo->DN(t,1);
 969     dt[0]=V1.X(); dt[1]=V1.Y();
 970   }
 971
 972   if(dtt){
 973     gp_Vec2d V2;
 974     V2=pargeo->DN(t,2);
 975     dtt[0]=V2.X(); dtt[1]=V2.Y();
 976   }
 977
 978   return 0;
 979 }
 980
 981 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 982                   real *duu, real *duv, real *dvv, void *user_data)
 983 {
 984   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
 985
 986   if(xyz){
 987    gp_Pnt P;
 988    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
 989    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
 990   }
 991
 992   if(du && dv){
 993     gp_Pnt P;
 994     gp_Vec D1U,D1V;
 995
 996     geometry->D1(uv[0],uv[1],P,D1U,D1V);
 997     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
 998     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
 999   }
1000
1001   if(duu && duv && dvv){
1002     gp_Pnt P;
1003     gp_Vec D1U,D1V;
1004     gp_Vec D2U,D2V,D2UV;
1005
1006     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
1007     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
1008     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
1009     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
1010   }
1011
1012   return 0;
1013 }
1014
1015
1016 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
1017 {
1018   if (face_id == 1) {
1019     if (my_u_min > uv[0]) {
1020       my_u_min = uv[0];
1021     }
1022     if (my_v_min > uv[1]) {
1023       my_v_min = uv[1];
1024     }
1025     if (my_u_max < uv[0]) {
1026       my_u_max = uv[0];
1027     }
1028     if (my_v_max < uv[1]) {
1029       my_v_max = uv[1];
1030     }
1031   }
1032
1033   if (FaceId2PythonSmp.count(face_id) != 0){
1034     PyObject * pyresult = NULL;
1035     PyObject* new_stderr = NULL;
1036     assert(Py_IsInitialized());
1037     PyGILState_STATE gstate;
1038     gstate = PyGILState_Ensure();
1039     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],"(f,f)",uv[0],uv[1]);
1040     double result;
1041     if ( pyresult == NULL){
1042       fflush(stderr);
1043       string err_description="";
1044       new_stderr = newPyStdOut(err_description);
1045       PySys_SetObject("stderr", new_stderr);
1046       PyErr_Print();
1047       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1048       Py_DECREF(new_stderr);
1049       MESSAGE("Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description);
1050       result = *((double*)user_data);
1051       }
1052     else {
1053       result = PyFloat_AsDouble(pyresult);
1054       Py_DECREF(pyresult);
1055     }
1056     *size = result;
1057     //MESSAGE("f(" << uv[0] << "," << uv[1] << ")" << " = " << result);
1058     PyGILState_Release(gstate);
1059   }
1060   else {
1061     *size = *((double*)user_data);
1062   }
1063   return STATUS_OK;
1064 }
1065
1066 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
1067 {
1068   if (EdgeId2PythonSmp.count(edge_id) != 0){
1069     PyObject * pyresult = NULL;
1070     PyObject* new_stderr = NULL;
1071     assert(Py_IsInitialized());
1072     PyGILState_STATE gstate;
1073     gstate = PyGILState_Ensure();
1074     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],"(f)",t);
1075     double result;
1076     if ( pyresult == NULL){
1077       fflush(stderr);
1078       string err_description="";
1079       new_stderr = newPyStdOut(err_description);
1080       PySys_SetObject("stderr", new_stderr);
1081       PyErr_Print();
1082       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1083       Py_DECREF(new_stderr);
1084       MESSAGE("Can't evaluate f(" << t << ")" << " error is " << err_description);
1085       result = *((double*)user_data);
1086       }
1087     else {
1088       result = PyFloat_AsDouble(pyresult);
1089       Py_DECREF(pyresult);
1090     }
1091     *size = result;
1092     PyGILState_Release(gstate);
1093   }
1094   else {
1095     *size = *((double*)user_data);
1096   }
1097   return STATUS_OK;
1098 }
1099
1100 status_t size_on_vertex(integer point_id, real *size, void *user_data)
1101 {
1102   if (VertexId2PythonSmp.count(point_id) != 0){
1103     PyObject * pyresult = NULL;
1104     PyObject* new_stderr = NULL;
1105     assert(Py_IsInitialized());
1106     PyGILState_STATE gstate;
1107     gstate = PyGILState_Ensure();
1108     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],"");
1109     double result;
1110     if ( pyresult == NULL){
1111       fflush(stderr);
1112       string err_description="";
1113       new_stderr = newPyStdOut(err_description);
1114       PySys_SetObject("stderr", new_stderr);
1115       PyErr_Print();
1116       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1117       Py_DECREF(new_stderr);
1118       MESSAGE("Can't evaluate f()" << " error is " << err_description);
1119       result = *((double*)user_data);
1120       }
1121     else {
1122       result = PyFloat_AsDouble(pyresult);
1123       Py_DECREF(pyresult);
1124     }
1125     *size = result;
1126     PyGILState_Release(gstate);
1127   }
1128   else {
1129     *size = *((double*)user_data);
1130   }
1131  return STATUS_OK;
1132 }
1133
1134 status_t message_callback(message_t *msg, void *user_data)
1135 {
1136   integer errnumber = 0;
1137   char *desc;
1138   message_get_number(msg, &errnumber);
1139   message_get_description(msg, &desc);
1140   if ( errnumber < 0 ) {
1141     string * error = (string*)user_data;
1142 //   if ( !error->empty() )
1143 //     *error += "\n";
1144     // remove ^A from the tail
1145     int len = strlen( desc );
1146     while (len > 0 && desc[len-1] != '\n')
1147       len--;
1148     error->append( desc, len );
1149   }
1150   else {
1151       std::cout << desc << std::endl;
1152   }
1153   return STATUS_OK;
1154 }
1155
1156
1157 //=============================================================================
1158 /*!
1159  *
1160  */
1161 //=============================================================================
1162 bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh& aMesh,
1163                                    const TopoDS_Shape& aShape,
1164                                    MapShapeNbElems& aResMap)
1165 {
1166   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
1167   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
1168   //int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
1169   //double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
1170   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
1171   if(_hypothesis) {
1172     _physicalMesh  = (int) _hypothesis->GetPhysicalMesh();
1173     _phySize       = _hypothesis->GetPhySize();
1174     //_geometricMesh = (int) hyp->GetGeometricMesh();
1175     //_angleMeshS    = hyp->GetAngleMeshS();
1176     _angleMeshC    = _hypothesis->GetAngleMeshC();
1177   }
1178
1179   bool IsQuadratic = false;
1180
1181   // ----------------
1182   // evaluate 1D
1183   // ----------------
1184   TopTools_DataMapOfShapeInteger EdgesMap;
1185   double fullLen = 0.0;
1186   double fullNbSeg = 0;
1187   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
1188     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1189     if( EdgesMap.IsBound(E) )
1190       continue;
1191     SMESH_subMesh *sm = aMesh.GetSubMesh(E);
1192     double aLen = SMESH_Algo::EdgeLength(E);
1193     fullLen += aLen;
1194     int nb1d = 0;
1195     if(_physicalMesh==1) {
1196        nb1d = (int)( aLen/_phySize + 1 );
1197     }
1198     else {
1199       // use geometry
1200       double f,l;
1201       Handle(Geom_Curve) C = BRep_Tool::Curve(E,f,l);
1202       double fullAng = 0.0;
1203       double dp = (l-f)/200;
1204       gp_Pnt P1,P2,P3;
1205       C->D0(f,P1);
1206       C->D0(f+dp,P2);
1207       gp_Vec V1(P1,P2);
1208       for(int j=2; j<=200; j++) {
1209         C->D0(f+dp*j,P3);
1210         gp_Vec V2(P2,P3);
1211         fullAng += fabs(V1.Angle(V2));
1212         V1 = V2;
1213         P2 = P3;
1214       }
1215       nb1d = (int)( fullAng/_angleMeshC + 1 );
1216     }
1217     fullNbSeg += nb1d;
1218     std::vector<int> aVec(17);
1219     for(int i=0; i<17; i++) aVec[i]=0;
1220     if( IsQuadratic > 0 ) {
1221       aVec[0] = 2*nb1d - 1;
1222       aVec[2] = nb1d;
1223     }
1224     else {
1225       aVec[0] = nb1d - 1;
1226       aVec[1] = nb1d;
1227     }
1228     aResMap.insert(std::make_pair(sm,aVec));
1229     EdgesMap.Bind(E,nb1d);
1230   }
1231   double ELen = fullLen/fullNbSeg;
1232   // ----------------
1233   // evaluate 2D
1234   // ----------------
1235   // try to evaluate as in MEFISTO
1236   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
1237     TopoDS_Face F = TopoDS::Face( exp.Current() );
1238     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1239     GProp_GProps G;
1240     BRepGProp::SurfaceProperties(F,G);
1241     double anArea = G.Mass();
1242     int nb1d = 0;
1243     for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) {
1244       nb1d += EdgesMap.Find(exp1.Current());
1245     }
1246     int nbFaces = (int) ( anArea/( ELen*ELen*sqrt(3.) / 4 ) );
1247     int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1248     std::vector<int> aVec(17);
1249     for(int i=0; i<17; i++) aVec[i]=0;
1250     if( IsQuadratic ) {
1251       int nb1d_in = (nbFaces*3 - nb1d) / 2;
1252       aVec[0] = nbNodes + nb1d_in;
1253       aVec[4] = nbFaces;
1254     }
1255     else {
1256       aVec[0] = nbNodes;
1257       aVec[3] = nbFaces;
1258     }
1259     aResMap.insert(std::make_pair(sm,aVec));
1260   }
1261
1262   // ----------------
1263   // evaluate 3D
1264   // ----------------
1265   GProp_GProps G;
1266   BRepGProp::VolumeProperties(aShape,G);
1267   double aVolume = G.Mass();
1268   double tetrVol = 0.1179*ELen*ELen*ELen;
1269   int nbVols = (int)aVolume/tetrVol;
1270   int nb1d_in = (int) ( ( nbVols*6 - fullNbSeg ) / 6 );
1271   std::vector<int> aVec(17);
1272   for(int i=0; i<17; i++) aVec[i]=0;
1273   if( IsQuadratic ) {
1274     aVec[0] = nb1d_in/3 + 1 + nb1d_in;
1275     aVec[9] = nbVols;
1276   }
1277   else {
1278     aVec[0] = nb1d_in/3 + 1;
1279     aVec[8] = nbVols;
1280   }
1281   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1282   aResMap.insert(std::make_pair(sm,aVec));
1283
1284   return true;
1285 }