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