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